Fitting the ST model for black carbon in Denver, CO

1 Introduction

In this script, we are building the spatiotemporal prediction model for black carbon in the Denver metro area. This model uses data from 5 sampling campaigns (as detailed in 15_Exploring_BC_Data.R). Two models are developed here: one which uses one temporal basis function (Model A) and one that uses two temporal basis functions (Model B).

Updates to this model are in response to reviewer comments on an earlier manuscript version submitted to ES&T.

A summary table at the end of the document summarizes the performance for each of these models. Overall, Model B performs well and allows us to predict BC concentrations at a weekly temporal resolution from 2009 through 2020.

## Loading required package: Matrix

2 Developing data sets

2.1 Denver data set

Here I’m reading in the filter data sets and getting them into the right shape for the SpatioTemporal package. Here the ST covariates cover the entire study period (2009 - 2020).

lur_data_all <- read_csv(here::here("Data", "Final_BC_Data_Set_for_ST_Model.csv"))

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   site_id = col_character(),
##   filter_id = col_character(),
##   campaign = col_character(),
##   StartDateTimeLocal = col_logical(),
##   EndDateTimeLocal = col_logical(),
##   sample_week = col_date(format = ""),
##   logged_runtime = col_logical(),
##   logged_rt_volume_L = col_logical(),
##   low_volume_flag = col_logical(),
##   ultralow_volume_flag = col_logical(),
##   pm_below_lod = col_logical(),
##   pm_below_loq = col_logical(),
##   bc_below_lod = col_logical(),
##   negative_pm_mass = col_logical(),
##   potential_contamination = col_logical(),
##   pm_mass_ug = col_logical(),
##   bc_mass_ug_corrected = col_logical(),
##   blank_corrected_bc = col_logical(),
##   bc_blank_mean = col_logical(),
##   pm_ug_m3_raw = col_logical()
##   # ... with 36 more columns
## )

## See spec(...) for full column specifications.

lur_data <- lur_data_all

bc_obs <- lur_data %>%
  dplyr::select(site_id, st_week, bc_ug_m3) %>%
  group_by(site_id, st_week) %>%
  summarize(bc_ug_m3 = mean(bc_ug_m3, na.rm = T)) %>%
  mutate(bc_ug_m3 = ifelse(is.nan(bc_ug_m3), NA, bc_ug_m3)) %>%
  mutate(log_bc = log(bc_ug_m3)) %>%
  dplyr::select(-bc_ug_m3) 

#' Pivot to a wide data frame
bc_obs2 <- bc_obs %>%
  pivot_wider(id_cols = st_week, names_from = site_id, values_from = log_bc) %>%
  # names_from = site_id, values_from = bc_ug_m3) %>%
  as.data.frame() %>%
  arrange(st_week)
rownames(bc_obs2) <- bc_obs2$st_week
bc_obs2$st_week <- NULL
bc_obs2 <- as.matrix(bc_obs2)

bc_weeks <- rownames(bc_obs2)
tail(bc_weeks)

## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"

#' Spatial covariates (scaled)
load(here::here("Results", "BC_LASSO_Model_5C.Rdata"))

lasso_coef_df <- data.frame(name = log_bc_lasso_coef3@Dimnames[[1]][log_bc_lasso_coef3@i + 1],
                            coefficient = log_bc_lasso_coef3@x)
covars <- as.character(lasso_coef_df$name[-1])
covars 

## [1] "impervious_2500" "open_2500"       "med_int_2500"    "high_int_100"   
## [5] "dist_m_cafos"    "aadt_100"        "aadt_250"

covar_fun <- paste("~", paste(covars, collapse = " + "))
covar_fun

## [1] "~ impervious_2500 + open_2500 + med_int_2500 + high_int_100 + dist_m_cafos + aadt_100 + aadt_250"

bc_sp_cov <- select(lur_data, site_id, lon, lat, all_of(covars)) %>%
  distinct() %>%
  mutate_at(.vars = vars(all_of(covars)),
            scale)
summary(bc_sp_cov)

##    site_id               lon              lat         impervious_2500.V1 
##  Length:69          Min.   :-105.2   Min.   :39.57   Min.   :-1.9791101  
##  Class :character   1st Qu.:-105.0   1st Qu.:39.72   1st Qu.:-0.7093639  
##  Mode  :character   Median :-104.9   Median :39.75   Median : 0.1278296  
##                     Mean   :-104.9   Mean   :39.76   Mean   : 0.0000000  
##                     3rd Qu.:-104.8   3rd Qu.:39.79   3rd Qu.: 0.4803932  
##                     Max.   :-104.7   Max.   :39.95   Max.   : 1.8291463  
##      open_2500.V1       med_int_2500.V1      high_int_100.V1  
##  Min.   :-1.3794337   Min.   :-2.1968789   Min.   :-0.532437  
##  1st Qu.:-0.5497931   1st Qu.:-0.7166467   1st Qu.:-0.532437  
##  Median :-0.1781892   Median : 0.2108832   Median :-0.532437  
##  Mean   : 0.0000000   Mean   : 0.0000000   Mean   : 0.000000  
##  3rd Qu.: 0.4072220   3rd Qu.: 0.7753690   3rd Qu.: 0.015894  
##  Max.   : 2.7693377   Max.   : 1.7719166   Max.   : 4.768092  
##    dist_m_cafos.V1        aadt_100.V1         aadt_250.V1    
##  Min.   :-2.1610758   Min.   :-0.390807   Min.   :-0.664113  
##  1st Qu.:-0.4079546   1st Qu.:-0.390807   1st Qu.:-0.664113  
##  Median :-0.0010458   Median :-0.390807   Median :-0.365247  
##  Mean   : 0.0000000   Mean   : 0.000000   Mean   : 0.000000  
##  3rd Qu.: 0.3418171   3rd Qu.:-0.115098   3rd Qu.: 0.063484  
##  Max.   : 2.2925064   Max.   : 4.632471   Max.   : 4.155107

bc_sp_cov <- bc_sp_cov %>%
  rename(ID = site_id) %>%
  mutate(lon2 = lon, lat2 = lat) %>%
  st_as_sf(coords = c("lon2", "lat2"), crs = ll_wgs84) %>%
  st_transform(crs = albers)
sp_coords <- do.call(rbind, st_geometry(bc_sp_cov)) %>% as_tibble()
names(sp_coords) <- c("x", "y")

bc_sp_cov <- bind_cols(bc_sp_cov, sp_coords) %>%
  st_set_geometry(NULL) %>%
  as.data.frame()
head(bc_sp_cov)

bc_sp_cov$type <- ifelse(bc_sp_cov$ID == "central", "central", "dist")
bc_sp_cov$type <- as.factor(bc_sp_cov$type)
summary(bc_sp_cov)

##       ID                 lon              lat         impervious_2500.V1 
##  Length:69          Min.   :-105.2   Min.   :39.57   Min.   :-1.9791101  
##  Class :character   1st Qu.:-105.0   1st Qu.:39.72   1st Qu.:-0.7093639  
##  Mode  :character   Median :-104.9   Median :39.75   Median : 0.1278296  
##                     Mean   :-104.9   Mean   :39.76   Mean   : 0.0000000  
##                     3rd Qu.:-104.8   3rd Qu.:39.79   3rd Qu.: 0.4803932  
##                     Max.   :-104.7   Max.   :39.95   Max.   : 1.8291463  
##      open_2500.V1       med_int_2500.V1      high_int_100.V1  
##  Min.   :-1.3794337   Min.   :-2.1968789   Min.   :-0.532437  
##  1st Qu.:-0.5497931   1st Qu.:-0.7166467   1st Qu.:-0.532437  
##  Median :-0.1781892   Median : 0.2108832   Median :-0.532437  
##  Mean   : 0.0000000   Mean   : 0.0000000   Mean   : 0.000000  
##  3rd Qu.: 0.4072220   3rd Qu.: 0.7753690   3rd Qu.: 0.015894  
##  Max.   : 2.7693377   Max.   : 1.7719166   Max.   : 4.768092  
##    dist_m_cafos.V1        aadt_100.V1         aadt_250.V1           x          
##  Min.   :-2.1610758   Min.   :-0.390807   Min.   :-0.664113   Min.   :-782614  
##  1st Qu.:-0.4079546   1st Qu.:-0.390807   1st Qu.:-0.664113   1st Qu.:-763240  
##  Median :-0.0010458   Median :-0.390807   Median :-0.365247   Median :-756412  
##  Mean   : 0.0000000   Mean   : 0.000000   Mean   : 0.000000   Mean   :-756776  
##  3rd Qu.: 0.3418171   3rd Qu.:-0.115098   3rd Qu.: 0.063484   3rd Qu.:-749434  
##  Max.   : 2.2925064   Max.   : 4.632471   Max.   : 4.155107   Max.   :-738740  
##        y                type   
##  Min.   :1873044   central: 1  
##  1st Qu.:1890776   dist   :68  
##  Median :1894956               
##  Mean   :1895217               
##  3rd Qu.:1899263               
##  Max.   :1917328

cor(bc_sp_cov[,covars])

##                 impervious_2500  open_2500 med_int_2500 high_int_100
## impervious_2500      1.00000000 -0.8409441    0.8987424   0.48006634
## open_2500           -0.84094408  1.0000000   -0.8002841  -0.37979903
## med_int_2500         0.89874242 -0.8002841    1.0000000   0.39325577
## high_int_100         0.48006634 -0.3797990    0.3932558   1.00000000
## dist_m_cafos        -0.02217726  0.3367603   -0.2416443  -0.01905821
## aadt_100             0.40911376 -0.3329512    0.2614633   0.71529092
## aadt_250             0.42742607 -0.3251319    0.2768334   0.60158329
##                 dist_m_cafos   aadt_100   aadt_250
## impervious_2500  -0.02217726  0.4091138  0.4274261
## open_2500         0.33676028 -0.3329512 -0.3251319
## med_int_2500     -0.24164427  0.2614633  0.2768334
## high_int_100     -0.01905821  0.7152909  0.6015833
## dist_m_cafos      1.00000000  0.1015568  0.2311201
## aadt_100          0.10155682  1.0000000  0.8062597
## aadt_250          0.23112012  0.8062597  1.0000000

#' NO2, and smoke spatiotemporal predictors
#' NO2 at each site estimated using IDW
bc_st_no2 <- dplyr::select(lur_data, site_id, st_week, idw_no2) %>%
  distinct() %>%
  arrange(st_week) %>%
  pivot_wider(id_cols = st_week,
              names_from = site_id, values_from = idw_no2) %>%
  as.data.frame()
rownames(bc_st_no2) <- bc_st_no2$st_week
bc_st_no2$st_week <- NULL
bc_st_no2 <- as.matrix(bc_st_no2)

no2_weeks <- rownames(bc_st_no2)
tail(no2_weeks)

## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"

setdiff(bc_weeks, no2_weeks)

## character(0)

#' Smoke day indicator variable based on WFS paper method (see Brey and Fischer, 2016)
#' based on any smoke variable in the area using the 2 sd cut-off (area_smoke_2sd)
bc_st_smk <- dplyr::select(lur_data, site_id, st_week, area_smoke_2sd) %>%
  distinct() %>%
  arrange(st_week) %>%
  pivot_wider(id_cols = st_week,
              names_from = site_id, values_from = area_smoke_2sd) %>%
  as.data.frame()
rownames(bc_st_smk) <- bc_st_smk$st_week
bc_st_smk$st_week <- NULL
bc_st_smk <- as.matrix(bc_st_smk)

smk_weeks <- rownames(bc_st_smk)
tail(smk_weeks)

## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"

tail(bc_weeks)

## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"

setdiff(bc_weeks, smk_weeks)

## character(0)

Create a new denver.data object.

denver.data <- createSTdata(obs = bc_obs2,
                            covars = bc_sp_cov,
                            SpatioTemporal = list(bc_st_no2 = bc_st_no2,
                                                  bc_st_smk = bc_st_smk))
D <- createDataMatrix(denver.data)
denver.data

## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 231 (observed: 231)
##  No. obs: 943
## 
## Only constant temporal trend,with dates:
##  2016-02-08 to 2020-11-30
## 
## 13 covariate(s):
##  [1] "ID"              "lon"             "lat"             "impervious_2500"
##  [5] "open_2500"       "med_int_2500"    "high_int_100"    "dist_m_cafos"   
##  [9] "aadt_100"        "aadt_250"        "x"               "y"              
## [13] "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 712
##   Dates: 2018-05-07 to 2020-04-06

2.2 Creating the NO2, BC, and PM data objects

Next we need to combine monitoring data for NO2, BC, and PM2.5 in order to calculate the time trends.

First up: NO2

Plot the time trends for NO2 at the monitors

Original units (ppb)

no2_ts <- ggplot(data = no2_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black", method = "gam") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
no2_ts

## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

Log-transformed NO2

no2_ts2 <- ggplot(data = no2_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black", method = "gam") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
no2_ts2

## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

Organize the observations and spatial data. Use the log-transformed NO2 here

#' NO2 observations
no2_obs <- no2_data %>%
  st_set_geometry(NULL) %>%
  dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
  mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
  filter(st_week %in% bc_weeks) %>%
  group_by(monitor_id, st_week) %>%
  mutate(log_mean = log(Arithmetic_Mean)) %>%
  summarize(no2 = mean(log_mean, na.rm=T)) %>%
  pivot_wider(id_cols = st_week,
              names_from = monitor_id, values_from = no2) %>%
  as.data.frame() %>%
  arrange(st_week)
head(no2_obs)

head(no2_obs$st_week)

## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"

tail(no2_obs$st_week)

## [1] "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23" "2020-11-30"
## [6] "2020-12-07"

# NO2 SP object
no2_data2 <- read_csv(here::here("Data", "Monitor_NO2_Data_AEA.csv")) %>%
  mutate(year = year(Date_Local)) %>%
  dplyr::select(-year) %>%
  st_as_sf(wkt = "WKT", crs = albers) %>%
  mutate(monitor_id = paste0(monitor_id, "_no2")) %>%
  filter(monitor_id %in% colnames(no2_obs))

## Parsed with column specification:
## cols(
##   .default = col_character(),
##   Parameter_Code = col_double(),
##   POC = col_double(),
##   Date_Local = col_date(format = ""),
##   Observation_Count = col_double(),
##   Observation_Percent = col_double(),
##   Arithmetic_Mean = col_double(),
##   Max_Value = col_double(),
##   `1st_Max_Hour` = col_double(),
##   AQI = col_double(),
##   Date_of_Last_Change = col_date(format = "")
## )

## See spec(...) for full column specifications.

no2_coords <- do.call(rbind, st_geometry(no2_data2)) %>% as_tibble()
names(no2_coords) <- c("x", "y")
no2_sp_lonlat <- no2_data2 %>%
  st_transform(crs = ll_wgs84)
no2_coords2 <- do.call(rbind, st_geometry(no2_sp_lonlat)) %>%
  as_tibble() %>% setNames(c("lon","lat"))

no2_sp_cov <- st_set_geometry(no2_data2, NULL) %>%
  dplyr::select(monitor_id) %>%
  rename(ID = monitor_id)
no2_sp_cov <- bind_cols(no2_sp_cov, no2_coords, no2_coords2) %>%
  as.data.frame() %>%
  distinct()

Next, BC at the central site monitor

#' BC central site concentrations
bc_cent_data <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
  filter(!is.na(Arithmetic_Mean)) %>%
  st_as_sf(wkt = "WKT", crs = albers) %>%
  mutate(County_Code = str_sub(monitor_id, start = 3, end = 5)) %>%
  filter(County_Code %in% counties) %>%
  arrange(Date_Local, monitor_id) %>%
  mutate(Arithmetic_Mean = ifelse(Arithmetic_Mean <= 0.005, NA, Arithmetic_Mean))

## Parsed with column specification:
## cols(
##   .default = col_character(),
##   Parameter_Code = col_double(),
##   POC = col_double(),
##   Pollutant_Standard = col_logical(),
##   Date_Local = col_date(format = ""),
##   Observation_Count = col_double(),
##   Observation_Percent = col_double(),
##   Arithmetic_Mean = col_double(),
##   Max_Value = col_double(),
##   `1st_Max_Hour` = col_double(),
##   AQI = col_logical(),
##   Method_Code = col_double(),
##   Date_of_Last_Change = col_date(format = "")
## )

## See spec(...) for full column specifications.

head(bc_cent_data$Date_Local)

## [1] "2016-02-08" "2016-02-09" "2016-02-10" "2016-02-11" "2016-02-12"
## [6] "2016-02-13"

tail(bc_cent_data$Date_Local)

## [1] "2020-11-25" "2020-11-26" "2020-11-27" "2020-11-28" "2020-11-29"
## [6] "2020-11-30"

#' Add an identifier to differentiate these measurements from collocated pollutants at the same site
bc_cent_data$monitor_id <- paste0(bc_cent_data$monitor_id, "_bc")
unique(bc_cent_data$monitor_id)

## [1] "080310027_bc"

Plot the time trends for BC at the central site monitor

Original units (ug/m3)

bc_ts <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black", method = "gam") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  #facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
bc_ts

## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

Log-transformed BC

bc_ts2 <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  #facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
bc_ts2

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Organize the observations and spatial data

Use log-transformed data here

#' BC observations
bc_cent_obs <- bc_cent_data %>%
  st_set_geometry(NULL) %>%
  dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
  mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
  filter(st_week %in% bc_weeks) %>%
  group_by(monitor_id, st_week) %>%
  mutate(log_mean = log(Arithmetic_Mean)) %>%
  summarize(bc = mean(log_mean, na.rm=T)) %>%
  pivot_wider(id_cols = st_week,
              names_from = monitor_id, values_from = bc) %>%
  as.data.frame() %>%
  arrange(st_week)
head(bc_cent_obs)

head(bc_cent_obs$st_week)

## [1] "2016-02-08" "2016-02-15" "2016-02-22" "2016-02-29" "2016-03-07"
## [6] "2016-03-14"

tail(bc_cent_obs$st_week)

## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"

# BC Central Site SP object
bc_cent_data2 <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
  mutate(year = year(Date_Local)) %>%
  dplyr::select(-year) %>%
  st_as_sf(wkt = "WKT", crs = albers) %>%
  mutate(monitor_id = paste0(monitor_id, "_bc")) %>%
  filter(monitor_id %in% colnames(bc_cent_obs))

## Parsed with column specification:
## cols(
##   .default = col_character(),
##   Parameter_Code = col_double(),
##   POC = col_double(),
##   Pollutant_Standard = col_logical(),
##   Date_Local = col_date(format = ""),
##   Observation_Count = col_double(),
##   Observation_Percent = col_double(),
##   Arithmetic_Mean = col_double(),
##   Max_Value = col_double(),
##   `1st_Max_Hour` = col_double(),
##   AQI = col_logical(),
##   Method_Code = col_double(),
##   Date_of_Last_Change = col_date(format = "")
## )

## See spec(...) for full column specifications.

bc_cent_coords <- do.call(rbind, st_geometry(bc_cent_data2)) %>% as_tibble()
names(bc_cent_coords) <- c("x", "y")
bc_cent_sp_lonlat <- bc_cent_data2 %>%
  st_transform(crs = ll_wgs84)
bc_cent_coords2 <- do.call(rbind, st_geometry(bc_cent_sp_lonlat)) %>%
  as_tibble() %>% setNames(c("lon","lat"))

bc_cent_sp_cov <- st_set_geometry(bc_cent_data2, NULL) %>%
  dplyr::select(monitor_id) %>%
  rename(ID = monitor_id)
bc_cent_sp_cov <- bind_cols(bc_cent_sp_cov, bc_cent_coords, bc_cent_coords2) %>%
  as.data.frame() %>%
  distinct()

Lastly, PM2.5 from select monitors with long-term records

Plot the time trends for PM2.5 at the monitors

Original units (ug/m3)

pm_ts <- ggplot(data = pm_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
pm_ts

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Log-transformed PM2.5

pm_ts2 <- ggplot(data = pm_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black", method = "gam") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
pm_ts2

## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

Organize the observations and spatial data

Use the log-transformed data here

#' PM2.5 observations
pm_obs <- pm_data %>%
  st_set_geometry(NULL) %>%
  dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
  mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
  filter(st_week %in% bc_weeks) %>%
  group_by(monitor_id, st_week) %>%
  mutate(log_mean = log(Arithmetic_Mean)) %>%
  summarize(pm = mean(log_mean, na.rm=T)) %>%
  pivot_wider(id_cols = st_week,
              names_from = monitor_id, values_from = pm) %>%
  as.data.frame() %>%
  arrange(st_week)
head(pm_obs)

head(pm_obs$st_week)

## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"

tail(pm_obs$st_week)

## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"

# PM SP object
pm_data2 <- read_csv(here::here("Data", "Monitor_PM_Data_AEA.csv")) %>%
  mutate(year = year(Date_Local)) %>%
  dplyr::select(-year) %>%
  st_as_sf(wkt = "WKT", crs = albers) %>%
  mutate(monitor_id = paste0(monitor_id, "_pm")) %>%
  filter(monitor_id %in% colnames(pm_obs))

## Parsed with column specification:
## cols(
##   .default = col_character(),
##   Parameter_Code = col_double(),
##   POC = col_double(),
##   Date_Local = col_date(format = ""),
##   Observation_Count = col_double(),
##   Observation_Percent = col_double(),
##   Arithmetic_Mean = col_double(),
##   Max_Value = col_double(),
##   `1st_Max_Hour` = col_double(),
##   AQI = col_double(),
##   Method_Code = col_double(),
##   Date_of_Last_Change = col_date(format = "")
## )

## See spec(...) for full column specifications.

pm_coords <- do.call(rbind, st_geometry(pm_data2)) %>% as_tibble()
names(pm_coords) <- c("x", "y")
pm_sp_lonlat <- pm_data2 %>%
  st_transform(crs = ll_wgs84)
pm_coords2 <- do.call(rbind, st_geometry(pm_sp_lonlat)) %>%
  as_tibble() %>% setNames(c("lon","lat"))

pm_sp_cov <- st_set_geometry(pm_data2, NULL) %>%
  dplyr::select(monitor_id) %>%
  rename(ID = monitor_id)
pm_sp_cov <- bind_cols(pm_sp_cov, pm_coords, pm_coords2) %>%
  as.data.frame() %>%
  distinct()

Combine the observations and the spatial covariates Scale the pollutant measurements

#' Join log-transformed observations
pol_obs <- data.frame(st_week = rownames(bc_obs2)) %>%
  left_join(no2_obs, by = "st_week") %>% 
  left_join(bc_cent_obs, by = "st_week") %>%
  left_join(pm_obs, by = "st_week")
glimpse(pol_obs)

## Rows: 627
## Columns: 18
## $ st_week         <chr> "2008-12-29", "2009-01-05", "2009-01-12", "2009-01-19…
## $ `080013001_no2` <dbl> 2.4313422, 2.3337153, 3.0953820, 2.4960668, 2.6445846…
## $ `080310002_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310026_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_bc`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080010006_pm`  <dbl> 1.394118, 1.247789, 1.934844, 2.298456, 1.598321, 1.7…
## $ `080010008_pm`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080050005_pm`  <dbl> 1.2490759, 0.6817687, 1.3697744, 1.7730816, 1.1939225…
## $ `080130003_pm`  <dbl> 0.9679299, 0.8243293, 2.1047287, 1.7212142, 1.9059908…
## $ `080130012_pm`  <dbl> 1.0919008, 0.8908546, 1.2511276, 1.2703657, 0.8047190…
## $ `080310002_pm`  <dbl> 1.335209, 1.188381, 1.767997, 2.027844, 1.486213, 1.7…
## $ `080310023_pm`  <dbl> 1.661197, 1.513911, 1.919262, 1.996528, 1.909773, 1.7…
## $ `080310025_pm`  <dbl> 1.1537863, 0.9359011, 1.3262687, 1.9969858, 1.6521301…
## $ `080310026_pm`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_pm`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_pm`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…

rownames(pol_obs) <- pol_obs$st_week
pol_obs$st_week <- NULL
pol_obs <- as.matrix(pol_obs)

#' Scale the measurements to avoid issues where units are different
#' Remember! These have been log-transformed before scaling
#' Scale the data by columns (monitors)
rnames <- rownames(pol_obs)
cnames <- colnames(pol_obs)

pol_obs <- apply(pol_obs, 2, scale)

rownames(pol_obs) <- rnames
colnames(pol_obs) <- cnames

head(rownames(pol_obs))

## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"

tail(rownames(pol_obs))

## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"

head(colnames(pol_obs))

## [1] "080013001_no2" "080310002_no2" "080310026_no2" "080310027_no2"
## [5] "080310028_no2" "080310027_bc"

class(pol_obs)

## [1] "matrix" "array"

dim(pol_obs)

## [1] 627  17

#' SP variables
pol_sp_cov <- bind_rows(no2_sp_cov, bc_cent_sp_cov, pm_sp_cov)
pol_sp_cov

Create the data object

pol_STdata <- createSTdata(obs = pol_obs,
                           covars = pol_sp_cov)
print(pol_STdata)

## STdata-object with:
##  No. locations: 17 (observed: 16)
##  No. time points: 624 (observed: 624)
##  No. obs: 6574
## 
## Only constant temporal trend,with dates:
##  2008-12-29 to 2020-12-07
## 
## 5 covariate(s):
## [1] "ID"  "x"   "y"   "lon" "lat"
## 
## No spatio-temporal covariates.

2.3 Generate the temporal trend

The cross validation results for generating the time trends suggest 2 basis functions is the way to go, but let’s see if that holds up with our data

D_pol <- createDataMatrix(pol_STdata)
#D_pol

n_years <- length(2009:2020)
n_years*4

## [1] 48

n_years*8

## [1] 96

#' Trying 4 df per year
pol.SVD.cv.4py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 4)
print(pol.SVD.cv.4py)

## Result of SVDsmoothCV, average of CV-statistics:
##                 MSE        R2          AIC         BIC
## n.basis.0 0.9975662 0.0000000    0.9987807    5.097602
## n.basis.1 0.8016075 0.1964669  -96.6267210  -88.429078
## n.basis.2 0.6359015 0.3625417 -225.5766119 -213.280148
## n.basis.3 0.6299681 0.3684906 -229.4640525 -213.068767
## n.basis.4 0.6089310 0.3896015 -240.3947203 -219.900613

plot(pol.SVD.cv.4py)

#' Trying 8 df per year
pol.SVD.cv.8py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 8)
print(pol.SVD.cv.8py)

## Result of SVDsmoothCV, average of CV-statistics:
##                 MSE        R2          AIC         BIC
## n.basis.0 0.9975662 0.0000000    0.9987807    5.097602
## n.basis.1 0.7615784 0.2365880 -121.3718836 -113.174241
## n.basis.2 0.5808908 0.4176809 -266.8631552 -254.566691
## n.basis.3 0.5724708 0.4261239 -272.6314965 -256.236211
## n.basis.4 0.5602312 0.4384013 -279.8920086 -259.397901

plot(pol.SVD.cv.8py)

2.4 Plot the time trends

I tried using 1 and 2 basis functions for the time trends. I also checked out 4 and 8 df per year

2.4.1 1 basis function

First, see what the time trends for the monitoring data look like with 1 basis function.

4 df per year

pol_STdata1 <- updateTrend(pol_STdata, n.basis = 1, df = n_years * 4)

## Replacing existing trend.

pol_STdata1

## STdata-object with:
##  No. locations: 17 (observed: 16)
##  No. time points: 624 (observed: 624)
##  No. obs: 6574
## 
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 5 covariate(s):
## [1] "ID"  "x"   "y"   "lon" "lat"
## 
## No spatio-temporal covariates.

head(pol_STdata1$trend)

tail(pol_STdata1$trend)

plot(pol_STdata1$trend$date, pol_STdata1$trend$V1, col = 1, pch = 16, cex = 0.5,
     xlab = "Date", ylab = "Pollutants (scaled)")
lines(pol_STdata1$trend$date, pol_STdata1$atrend$V1, col = 1)

8 df per year

pol_STdata1a <- updateTrend(pol_STdata, n.basis = 1, df = n_years * 8)

## Replacing existing trend.

pol_STdata1a

## STdata-object with:
##  No. locations: 17 (observed: 16)
##  No. time points: 624 (observed: 624)
##  No. obs: 6574
## 
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 5 covariate(s):
## [1] "ID"  "x"   "y"   "lon" "lat"
## 
## No spatio-temporal covariates.

head(pol_STdata1a$trend)

tail(pol_STdata1a$trend)

plot(pol_STdata1a$trend$date, pol_STdata1a$trend$V1, col = 1, pch = 16, cex = 0.5,
     xlab = "Date", ylab = "Pollutants (scaled)")
lines(pol_STdata1a$trend$date, pol_STdata1a$atrend$V1, col = 1)

2.4.2 2 basis functions

Next, let’s see what 2 basis functions look like:

4 df per year

pol_STdata2 <- updateTrend(pol_STdata, n.basis = 2, df = n_years * 4)

## Replacing existing trend.

pol_STdata2

## STdata-object with:
##  No. locations: 17 (observed: 16)
##  No. time points: 624 (observed: 624)
##  No. obs: 6574
## 
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 5 covariate(s):
## [1] "ID"  "x"   "y"   "lon" "lat"
## 
## No spatio-temporal covariates.

head(pol_STdata2$trend)

tail(pol_STdata2$trend)

plot(pol_STdata2$trend$date, pol_STdata2$trend$V1, col = 1, pch = 16, cex = 0.5,
     xlab = "Date", ylab = "Pollutants (scaled)")
points(pol_STdata2$trend$date, pol_STdata2$trend$V2, col = 2, pch = 16, cex = 0.5)
lines(pol_STdata2$trend$date, pol_STdata2$trend$V1, col = 1)
lines(pol_STdata2$trend$date, pol_STdata2$trend$V2, col = 2)

8 df per year

pol_STdata2a <- updateTrend(pol_STdata, n.basis = 2, df = n_years * 8)

## Replacing existing trend.

pol_STdata2a

## STdata-object with:
##  No. locations: 17 (observed: 16)
##  No. time points: 624 (observed: 624)
##  No. obs: 6574
## 
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 5 covariate(s):
## [1] "ID"  "x"   "y"   "lon" "lat"
## 
## No spatio-temporal covariates.

head(pol_STdata2a$trend)

tail(pol_STdata2a$trend)

plot(pol_STdata2a$trend$date, pol_STdata2a$trend$V1, col = 1, pch = 16, cex = 0.5,
     xlab = "Date", ylab = "Pollutants (scaled)")
points(pol_STdata2a$trend$date, pol_STdata2a$trend$V2, col = 2, pch = 16, cex = 0.5)
lines(pol_STdata2a$trend$date, pol_STdata2a$trend$V1, col = 1)
lines(pol_STdata2a$trend$date, pol_STdata2a$trend$V2, col = 2)

2.5 Update the Denver data set

The next step is to update the trend data for the denver.data object and see if these trend data fit our observations.

Comparing the plots using one basis function vs. two basis functions, it looks like using two basis functions might be better than 1. There was some residual autocorrelation in the central site monitor data when using just one basis function (see the fourth plot in the first series of plots below). Using 8 degrees of freedom did not help.

Using two basis functions helped somewhat to reduce this autocorrelation, but it still persists. There are still some noticable trends in the residuals for the central monitoring site

2.5.1 1 basis function, 4 df per year

First update the trend data, then plot the trends for three distributed sites and the central BC monitoring site.

denver.data2.1 <- denver.data
denver.data2.1$trend <- pol_STdata1$trend
print(denver.data2.1)

## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 624 (observed: 231)
##  No. obs: 943
## 
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 13 covariate(s):
##  [1] "ID"              "lon"             "lat"             "impervious_2500"
##  [5] "open_2500"       "med_int_2500"    "high_int_100"    "dist_m_cafos"   
##  [9] "aadt_100"        "aadt_250"        "x"               "y"              
## [13] "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 712
##   Dates: 2018-05-07 to 2020-04-06

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_1")
plot(denver.data2.1, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_1")
plot(denver.data2.1, "pacf", ID="d_1")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_10")
plot(denver.data2.1, "pacf", ID="d_10")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_43")
plot(denver.data2.1, "pacf", ID="d_43")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_53")
plot(denver.data2.1, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_53")
plot(denver.data2.1, "pacf", ID="d_53")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend central")
plot(denver.data2.1, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="central")
plot(denver.data2.1, "pacf", ID="central")

2.5.2 1 basis function, 8 df per year

denver.data2.1a <- denver.data
denver.data2.1a$trend <- pol_STdata1a$trend
print(denver.data2.1a)

## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 624 (observed: 231)
##  No. obs: 943
## 
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 13 covariate(s):
##  [1] "ID"              "lon"             "lat"             "impervious_2500"
##  [5] "open_2500"       "med_int_2500"    "high_int_100"    "dist_m_cafos"   
##  [9] "aadt_100"        "aadt_250"        "x"               "y"              
## [13] "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 712
##   Dates: 2018-05-07 to 2020-04-06

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_1")
plot(denver.data2.1a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_1")
plot(denver.data2.1a, "pacf", ID="d_1")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_10")
plot(denver.data2.1a, "pacf", ID="d_10")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_43")
plot(denver.data2.1a, "pacf", ID="d_43")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_53")
plot(denver.data2.1a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_53")
plot(denver.data2.1a, "pacf", ID="d_53")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend central")
plot(denver.data2.1a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="central")
plot(denver.data2.1a, "pacf", ID="central")

2.5.3 2 basis functions, 4 df per year

denver.data2.2 <- denver.data
denver.data2.2$trend <- pol_STdata2$trend
print(denver.data2.2)

## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 624 (observed: 231)
##  No. obs: 943
## 
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 13 covariate(s):
##  [1] "ID"              "lon"             "lat"             "impervious_2500"
##  [5] "open_2500"       "med_int_2500"    "high_int_100"    "dist_m_cafos"   
##  [9] "aadt_100"        "aadt_250"        "x"               "y"              
## [13] "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 712
##   Dates: 2018-05-07 to 2020-04-06

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_1")
plot(denver.data2.2, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_1")
plot(denver.data2.2, "pacf", ID="d_1")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_10")
plot(denver.data2.2, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_10")
plot(denver.data2.2, "pacf", ID="d_10")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.2, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_43")
plot(denver.data2.2, "pacf", ID="d_43")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_53")
plot(denver.data2.2, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_53")
plot(denver.data2.2, "pacf", ID="d_53")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend central")
plot(denver.data2.2, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="central")
plot(denver.data2.2, "pacf", ID="central")

2.5.4 2 basis functions, 8 df per year

denver.data2.2a <- denver.data
denver.data2.2a$trend <- pol_STdata2a$trend
print(denver.data2.2a)

## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 624 (observed: 231)
##  No. obs: 943
## 
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 13 covariate(s):
##  [1] "ID"              "lon"             "lat"             "impervious_2500"
##  [5] "open_2500"       "med_int_2500"    "high_int_100"    "dist_m_cafos"   
##  [9] "aadt_100"        "aadt_250"        "x"               "y"              
## [13] "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 712
##   Dates: 2018-05-07 to 2020-04-06

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_1")
plot(denver.data2.2a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_1")
plot(denver.data2.2a, "pacf", ID="d_1")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_10")
plot(denver.data2.2a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_10")
plot(denver.data2.2a, "pacf", ID="d_10")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.2a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_43")
plot(denver.data2.2a, "pacf", ID="d_43")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_53")
plot(denver.data2.2a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_53")
plot(denver.data2.2a, "pacf", ID="d_53")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend central")
plot(denver.data2.2a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="central")
plot(denver.data2.2a, "pacf", ID="central")

3 Testing possible models

Building on the results of the previous modeling efforts, I’m going to update “Model 4.2”, which used an ’exp covariance structure for the error term and iid for the beta fields.

3.1 Model A: 1 temporal trend (basis) function (4 df per year)

3.1.1 Create the model object

For this version of the model, use iid for cov.beta (beta, beta1) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.

denver.data.A <- denver.data2.1
names(denver.data.A$covars)

##  [1] "ID"              "lon"             "lat"             "impervious_2500"
##  [5] "open_2500"       "med_int_2500"    "high_int_100"    "dist_m_cafos"   
##  [9] "aadt_100"        "aadt_250"        "x"               "y"              
## [13] "type"

LUR.A <- list(covar_fun, covar_fun)

cov.beta.A <-  list(covf = c("iid", "iid"), nugget = T)
cov.nu.A <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.A <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")

denver.model.A <- createSTmodel(denver.data.A, LUR = LUR.A,
                                ST = c("bc_st_no2", "bc_st_smk"),
                                cov.beta = cov.beta.A, cov.nu = cov.nu.A,
                                locations = locations.A)

## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.

denver.model.A

## STmodel-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 624 (observed: 231)
##  No. obs: 943
## 
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## Models for the beta-fields are:
## $const
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 + 
##     dist_m_cafos + aadt_100 + aadt_250
## 
## $V1
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 + 
##     dist_m_cafos + aadt_100 + aadt_250
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## Covariance model for the beta-field(s):
##  Covariance type(s): iid, iid 
##  Nugget: Yes, Yes 
## Covariance model for the nu-field(s):
##  Covariance type: exp 
##  Nugget: ~1 
##  Random effect: No 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 712
##   Dates: 2018-05-07 to 2020-04-06

3.1.2 Estimate model parameters

• nugget values can range from -1 to -6
• range values can range from 0 to 4
• try starting values where the sill and nugget are the same (e.g., both 0) and where they are very different (e.g., 0 and -5)
• make sure the starting values are pretty different

set.seed(1000)

names <- loglikeSTnames(denver.model.A, all=FALSE)
names

## [1] "log.nugget.const.iid"          "log.nugget.V1.iid"            
## [3] "nu.log.range.exp"              "nu.log.sill.exp"              
## [5] "nu.log.nugget.(Intercept).exp"

x.init.A <- cbind(c(0, 0, 0, 0, 0),
                  c(-1, -1, 4, -5, -1),
                  c(-5, -5, 4, -1, -5),
                  c(-1, -1, 8, -5, -1),
                  c(-5, -5, 8, -1, -5),
                  c(-1, -1, 12, -5, -1),
                  c(-5, -5, 12, -1, -5),
                  c(-7, -7, 12, -3, -3))

rownames(x.init.A) <- loglikeSTnames(denver.model.A, all=FALSE)
x.init.A

##                               [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## log.nugget.const.iid             0   -1   -5   -1   -5   -1   -5   -7
## log.nugget.V1.iid                0   -1   -5   -1   -5   -1   -5   -7
## nu.log.range.exp                 0    4    4    8    8   12   12   12
## nu.log.sill.exp                  0   -5   -1   -5   -1   -5   -1   -3
## nu.log.nugget.(Intercept).exp    0   -1   -5   -1   -5   -1   -5   -3

#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.A <- estimate.STmodel(denver.model.A, x.init.A)

## Optimisation using starting value 1/8

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=       428.12  |proj g|=           15
## At iterate    10  f =      -1187.6  |proj g|=        1.5251
## 
## iterations 18
## function evaluations 36
## segments explored during Cauchy searches 22
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0128187
## final function value -1187.79
## 
## F = -1187.79
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1187.786129 
## converged

## Optimisation using starting value 2/8

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=       -299.6  |proj g|=           14
## At iterate    10  f =      -1184.3  |proj g|=         0.423
## ys=-1.276e-02  -gs= 2.763e-02, BFGS update SKIPPED
## At iterate    20  f =      -1184.9  |proj g|=        11.462
## At iterate    30  f =      -1187.8  |proj g|=       0.19815
## ys=-8.732e-01  -gs= 1.086e+00, BFGS update SKIPPED
## At iterate    40  f =      -1403.6  |proj g|=       0.18215
## 
## iterations 43
## function evaluations 76
## segments explored during Cauchy searches 46
## BFGS updates skipped 2
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00508496
## final function value -1403.62
## 
## F = -1403.62
## final  value -1403.623292 
## converged

## Optimisation using starting value 3/8

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -421.29  |proj g|=           14
## At iterate    10  f =      -1183.3  |proj g|=        1.6336
## At iterate    20  f =      -1187.7  |proj g|=        0.6575
## ys=-7.955e+01  -gs= 7.565e-01, BFGS update SKIPPED
## At iterate    30  f =        -1401  |proj g|=          10.6
## At iterate    40  f =      -1409.4  |proj g|=       0.36759
## 
## iterations 48
## function evaluations 87
## segments explored during Cauchy searches 49
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0252493
## final function value -1409.4
## 
## F = -1409.4
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1409.403118 
## converged

## Optimisation using starting value 4/8

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -299.84  |proj g|=           14
## At iterate    10  f =      -1400.2  |proj g|=         7.974
## At iterate    20  f =        -1405  |proj g|=        7.8782
## At iterate    30  f =      -1407.6  |proj g|=        1.8038
## At iterate    40  f =      -1409.3  |proj g|=        2.6564
## At iterate    50  f =      -1409.4  |proj g|=    0.00075243
## 
## iterations 50
## function evaluations 67
## segments explored during Cauchy searches 53
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00075243
## final function value -1409.4
## 
## F = -1409.4
## final  value -1409.403111 
## converged

## Optimisation using starting value 5/8

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -495.77  |proj g|=           14
## At iterate    10  f =      -1409.3  |proj g|=       0.80461
## At iterate    20  f =      -1409.4  |proj g|=       0.15163
## 
## iterations 24
## function evaluations 36
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0995913
## final function value -1409.4
## 
## F = -1409.4
## final  value -1409.403060 
## converged

## Optimisation using starting value 6/8

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -302.63  |proj g|=           14
## At iterate    10  f =      -1408.9  |proj g|=        8.1359
## At iterate    20  f =      -1409.4  |proj g|=     0.0029418
## 
## iterations 21
## function evaluations 29
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.010386
## final function value -1409.4
## 
## F = -1409.4
## final  value -1409.403107 
## converged

## Optimisation using starting value 7/8

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1152.5  |proj g|=           14
## At iterate    10  f =      -1409.4  |proj g|=       0.33301
## At iterate    20  f =      -1409.4  |proj g|=      0.025957
## 
## iterations 20
## function evaluations 40
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.025957
## final function value -1409.4
## 
## F = -1409.4
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1409.403103 
## converged

## Optimisation using starting value 8/8

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1118.9  |proj g|=           12
## At iterate    10  f =      -1409.4  |proj g|=       0.13337
## 
## iterations 15
## function evaluations 23
## segments explored during Cauchy searches 18
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.127558
## final function value -1409.4
## 
## F = -1409.4
## final  value -1409.402987 
## converged

print(est.denver.model.A)

## Optimisation for STmodel with 8 starting points.
##   Results: 6 converged, 2 not converged, 0 failed.
##   Best result for starting point 3, optimisation has converged
## 
## No fixed parameters.
## 
## Estimated parameters for all starting point(s):
##                                         [,1]           [,2]         [,3]
## gamma.bc_st_no2                 0.0190177471   0.0139611091  0.013524343
## gamma.bc_st_smk                -0.0426958408  -0.0160214985 -0.017276051
## alpha.const.(Intercept)        -0.1503355265   0.0025614497  0.008365957
## alpha.const.impervious_2500     0.0376285789   0.0297836107  0.030578606
## alpha.const.open_2500          -0.0076291505  -0.0104722683 -0.010608880
## alpha.const.med_int_2500       -0.0164412434  -0.0205190931 -0.020872498
## alpha.const.high_int_100        0.0113010906   0.0123640421  0.010811603
## alpha.const.dist_m_cafos       -0.0373287719  -0.0354396137 -0.035629090
## alpha.const.aadt_100            0.0322796488   0.0292214710  0.029130554
## alpha.const.aadt_250            0.0164520345   0.0148726789  0.014634801
## alpha.V1.(Intercept)            0.1874787211   0.2202093765  0.225017101
## alpha.V1.impervious_2500        0.0117787671  -0.0089070999 -0.009956409
## alpha.V1.open_2500              0.0118058417   0.0114616390  0.011181075
## alpha.V1.med_int_2500           0.0090074275   0.0130834301  0.011830049
## alpha.V1.high_int_100           0.0180754356   0.0178398056  0.018430627
## alpha.V1.dist_m_cafos          -0.0002816461  -0.0004239591 -0.002006719
## alpha.V1.aadt_100              -0.0323430204  -0.0319474609 -0.030466145
## alpha.V1.aadt_250              -0.0150887087  -0.0127119439 -0.010031267
## log.nugget.const.iid          -14.7001121803 -14.9989138692 -7.558680048
## log.nugget.V1.iid             -15.0000000000 -14.9994643477 -7.132387053
## nu.log.range.exp                0.0000000000  12.5059102289 12.520026954
## nu.log.sill.exp                -4.3239982804  -3.0584488703 -3.032854210
## nu.log.nugget.(Intercept).exp  -4.1033971401  -4.6498599256 -4.759856021
##                                       [,4]         [,5]         [,6]
## gamma.bc_st_no2                0.013525741  0.013529300  0.013525486
## gamma.bc_st_smk               -0.017271850 -0.017263869 -0.017272711
## alpha.const.(Intercept)        0.008340578  0.008261481  0.008346703
## alpha.const.impervious_2500    0.030578311  0.030574416  0.030578433
## alpha.const.open_2500         -0.010609062 -0.010610972 -0.010609050
## alpha.const.med_int_2500      -0.020872571 -0.020871243 -0.020872644
## alpha.const.high_int_100       0.010812892  0.010811989  0.010812616
## alpha.const.dist_m_cafos      -0.035630328 -0.035631806 -0.035630314
## alpha.const.aadt_100           0.029131348  0.029132930  0.029131662
## alpha.const.aadt_250           0.014630979  0.014625332  0.014630607
## alpha.V1.(Intercept)           0.225008013  0.224984562  0.225008189
## alpha.V1.impervious_2500      -0.009955584 -0.009961123 -0.009956377
## alpha.V1.open_2500             0.011183188  0.011185787  0.011183151
## alpha.V1.med_int_2500          0.011827663  0.011826552  0.011827449
## alpha.V1.high_int_100          0.018433956  0.018435501  0.018434386
## alpha.V1.dist_m_cafos         -0.002006304 -0.002010938 -0.002006642
## alpha.V1.aadt_100             -0.030462478 -0.030452721 -0.030461775
## alpha.V1.aadt_250             -0.010032886 -0.010031414 -0.010032386
## log.nugget.const.iid          -7.560464968 -7.560107182 -7.560275086
## log.nugget.V1.iid             -7.131144030 -7.127587858 -7.130684472
## nu.log.range.exp              12.519691760 12.517936372 12.519742836
## nu.log.sill.exp               -3.033134958 -3.033381695 -3.033081239
## nu.log.nugget.(Intercept).exp -4.759862263 -4.759736916 -4.759927367
##                                       [,7]         [,8]
## gamma.bc_st_no2                0.013523845  0.013523998
## gamma.bc_st_smk               -0.017273875 -0.017269878
## alpha.const.(Intercept)        0.008380116  0.008390079
## alpha.const.impervious_2500    0.030583778  0.030585150
## alpha.const.open_2500         -0.010607883 -0.010605271
## alpha.const.med_int_2500      -0.020873883 -0.020875042
## alpha.const.high_int_100       0.010808899  0.010821491
## alpha.const.dist_m_cafos      -0.035630538 -0.035630373
## alpha.const.aadt_100           0.029132549  0.029126451
## alpha.const.aadt_250           0.014630832  0.014637730
## alpha.V1.(Intercept)           0.225021583  0.225035682
## alpha.V1.impervious_2500      -0.009955880 -0.009932152
## alpha.V1.open_2500             0.011181831  0.011183267
## alpha.V1.med_int_2500          0.011824019  0.011826210
## alpha.V1.high_int_100          0.018438465  0.018434686
## alpha.V1.dist_m_cafos         -0.002005353 -0.001990822
## alpha.V1.aadt_100             -0.030465340 -0.030482658
## alpha.V1.aadt_250             -0.010026484 -0.010046108
## log.nugget.const.iid          -7.557186151 -7.567336821
## log.nugget.V1.iid             -7.127421757 -7.140153420
## nu.log.range.exp              12.520867210 12.522302555
## nu.log.sill.exp               -3.032945560 -3.033012579
## nu.log.nugget.(Intercept).exp -4.759776099 -4.758603669
## 
## Function value(s):
## [1] 1187.786 1403.623 1409.403 1409.403 1409.403 1409.403 1409.403 1409.403

3.1.3 Cross-validation

Define the CV groups

set.seed(1000)
site_idsA <- unique(denver.model.A$obs$ID[which(str_detect(denver.model.A$obs$ID, "d_") == T)])
site_idsA

##  [1] "d_1"  "d_18" "d_22" "d_25" "d_28" "d_31" "d_33" "d_36" "d_39" "d_4" 
## [11] "d_7"  "d_11" "d_12" "d_19" "d_23" "d_26" "d_29" "d_32" "d_35" "d_38"
## [21] "d_40" "d_8"  "d_42" "d_44" "d_46" "d_50" "d_52" "d_54" "d_56" "d_58"
## [31] "d_61" "d_63" "d_65" "d_68" "d_15" "d_41" "d_43" "d_45" "d_47" "d_51"
## [41] "d_53" "d_67" "d_49" "d_14" "d_48" "d_37" "d_55" "d_57" "d_64" "d_66"
## [51] "d_16" "d_2"  "d_5"  "d_60" "d_62" "d_34" "d_10" "d_13" "d_21" "d_3" 
## [61] "d_6"  "d_17" "d_20"

Ind.cv.A <- createCV(denver.model.A, groups = 10, min.dist = .1,
                     subset = site_idsA)

ID.cv.A <- sapply(split(denver.model.A$obs$ID, Ind.cv.A), unique)
print(sapply(ID.cv.A, length))

##  0  1  2  3  4  5  6  7  8  9 10 
##  1  7  7  7  6  6  6  6  6  6  6

table(Ind.cv.A)

## Ind.cv.A
##   0   1   2   3   4   5   6   7   8   9  10 
## 231  66  76  82  67  57  85  72  78  54  75

I.col.A <- apply(sapply(ID.cv.A, function(x) denver.model.A$locations$ID%in% x), 1,
                       function(x) if(sum(x)==1) which(x) else 0)
names(I.col.A) <- denver.model.A$locations$ID
print(I.col.A)

## central     d_1    d_10    d_11    d_12    d_13    d_14    d_15    d_16    d_17 
##       1       6       2       3       4       4       4       9       4      10 
##    d_18    d_19     d_2    d_20    d_21    d_22    d_23    d_25    d_26    d_28 
##       8       5      11       3       2      11       4       9      11       9 
##    d_29     d_3    d_31    d_32    d_33    d_34    d_35    d_36    d_37    d_38 
##       2       4       6       8       5       5      10       8       3       8 
##    d_39     d_4    d_40    d_41    d_42    d_43    d_44    d_45    d_46    d_47 
##       6       7       7       3       7       8      11       3       9       5 
##    d_48    d_49     d_5    d_50    d_51    d_52    d_53    d_54    d_55    d_56 
##       9       4       6       9      10       3       7       5       7       6 
##    d_57    d_58     d_6    d_60    d_61    d_62    d_63    d_64    d_65    d_66 
##       2       2      11       5       2      11      10       3       8       6 
##    d_67    d_68     d_7     d_8    d_24    d_27    d_30    d_59     d_9 
##       7      10       2      10       0       0       0       0       0

par(mfrow=c(1,1))
plot(denver.model.A$locations$long,
     denver.model.A$locations$lat,
     pch=23+floor(I.col.A/max(I.col.A)+.5), bg=I.col.A,
     xlab="Longitude", ylab="Latitude")

map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)

## [[1]]
## NULL

ID starting conditions using previous model without CV:

#' ID starting conditions using previous model without CV:
x.init.A.cv <- coef(est.denver.model.A, pars="cov")[,c("par","init")]
x.init.A.cv

Run the model with cross validation

est.denver.A.cv <- estimateCV(denver.model.A, x.init.A.cv, Ind.cv.A)

## 
## ***************************

## Estimation of CV-set 1/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1279.4  |proj g|=       13.463
## At iterate    10  f =      -1279.9  |proj g|=       0.17604
## 
## iterations 16
## function evaluations 24
## segments explored during Cauchy searches 16
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00869475
## final function value -1279.92
## 
## F = -1279.92
## final  value -1279.916549 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -391.08  |proj g|=           14
## At iterate    10  f =        -1086  |proj g|=        2.3228
## ys=-4.014e-01  -gs= 1.977e-01, BFGS update SKIPPED
## At iterate    20  f =      -1258.2  |proj g|=        19.721
## At iterate    30  f =      -1275.4  |proj g|=     0.0050604
## 
## iterations 30
## function evaluations 49
## segments explored during Cauchy searches 31
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00506041
## final function value -1275.42
## 
## F = -1275.42
## final  value -1275.424691 
## converged

## 

## 
## ***************************

## Estimation of CV-set 2/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1280.8  |proj g|=       1.2959
## At iterate    10  f =      -1280.9  |proj g|=       0.40948
## 
## iterations 15
## function evaluations 31
## segments explored during Cauchy searches 15
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0212405
## final function value -1280.93
## 
## F = -1280.93
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1280.932738 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -387.84  |proj g|=           14
## At iterate    10  f =      -1086.5  |proj g|=        2.4518
## At iterate    20  f =      -1090.4  |proj g|=        3.3494
## At iterate    30  f =      -1272.2  |proj g|=        3.3665
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1275.520364 
## stopped after 37 iterations

## 

## 
## ***************************

## Estimation of CV-set 3/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1285.1  |proj g|=        10.24
## 
## iterations 9
## function evaluations 18
## segments explored during Cauchy searches 11
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0629467
## final function value -1285.94
## 
## F = -1285.94
## final  value -1285.938947 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -378.99  |proj g|=           14
## At iterate    10  f =      -1080.4  |proj g|=         13.47
## ys=-3.594e-01  -gs= 6.181e-01, BFGS update SKIPPED
## ys=-2.044e+00  -gs= 2.394e+00, BFGS update SKIPPED
## At iterate    20  f =        -1196  |proj g|=        18.249
## At iterate    30  f =      -1276.8  |proj g|=        12.846
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1278.544393 
## stopped after 37 iterations

## 

## 
## ***************************

## Estimation of CV-set 4/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1273.2  |proj g|=       16.881
## At iterate    10  f =        -1274  |proj g|=       0.15717
## 
## iterations 17
## function evaluations 34
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0382929
## final function value -1274.04
## 
## F = -1274.04
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1274.042573 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -390.46  |proj g|=           14
## At iterate    10  f =        -1083  |proj g|=        2.3331
## At iterate    20  f =      -1087.3  |proj g|=       0.31284
## ys=-6.049e-01  -gs= 1.042e+00, BFGS update SKIPPED
## At iterate    30  f =      -1147.2  |proj g|=          20.9
## At iterate    40  f =        -1269  |proj g|=        10.561
## 
## iterations 45
## function evaluations 65
## segments explored during Cauchy searches 46
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00406416
## final function value -1270.08
## 
## F = -1270.08
## final  value -1270.080154 
## converged

## 

## 
## ***************************

## Estimation of CV-set 5/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1300.8  |proj g|=       9.1053
## At iterate    10  f =      -1301.2  |proj g|=      0.080753
## 
## iterations 17
## function evaluations 40
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0469096
## final function value -1301.24
## 
## F = -1301.24
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1301.239089 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -395.42  |proj g|=           14
## At iterate    10  f =      -1100.6  |proj g|=        1.8267
## At iterate    20  f =        -1105  |proj g|=       0.36357
## ys=-2.770e+00  -gs= 4.355e-02, BFGS update SKIPPED
## At iterate    30  f =      -1292.1  |proj g|=        19.885
## 
## iterations 38
## function evaluations 50
## segments explored during Cauchy searches 39
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0040607
## final function value -1297.73
## 
## F = -1297.73
## final  value -1297.728676 
## converged

## 

## 
## ***************************

## Estimation of CV-set 6/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1273.3  |proj g|=       8.8806
## At iterate    10  f =      -1273.6  |proj g|=      0.035777
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## 
## iterations 16
## function evaluations 42
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00197049
## final function value -1273.62
## 
## F = -1273.62
## final  value -1273.617044 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -382.95  |proj g|=           14
## At iterate    10  f =        -1065  |proj g|=        3.3156
## At iterate    20  f =      -1069.2  |proj g|=       0.27056
## At iterate    30  f =      -1267.7  |proj g|=        5.0383
## 
## iterations 36
## function evaluations 61
## segments explored during Cauchy searches 37
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0236423
## final function value -1268.61
## 
## F = -1268.61
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1268.613291 
## converged

## 

## 
## ***************************

## Estimation of CV-set 7/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1281.2  |proj g|=        7.644
## At iterate    10  f =      -1282.8  |proj g|=       0.13879
## At iterate    20  f =      -1282.8  |proj g|=      0.039849
## 
## iterations 20
## function evaluations 37
## segments explored during Cauchy searches 20
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0398489
## final function value -1282.79
## 
## F = -1282.79
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1282.788383 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -388.62  |proj g|=           14
## At iterate    10  f =      -1080.1  |proj g|=        2.3357
## ys=-4.483e-02  -gs= 2.380e-01, BFGS update SKIPPED
## ys=-1.945e+02  -gs= 6.398e+00, BFGS update SKIPPED
## At iterate    20  f =      -1260.6  |proj g|=        8.8126
## At iterate    30  f =        -1280  |proj g|=       0.19367
## 
## iterations 32
## function evaluations 43
## segments explored during Cauchy searches 33
## BFGS updates skipped 2
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.116086
## final function value -1280.01
## 
## F = -1280.01
## final  value -1280.007071 
## converged

## 

## 
## ***************************

## Estimation of CV-set 8/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1263.8  |proj g|=       9.6277
## At iterate    10  f =        -1265  |proj g|=       0.18815
## At iterate    20  f =        -1265  |proj g|=      0.015515
## 
## iterations 22
## function evaluations 37
## segments explored during Cauchy searches 22
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0126596
## final function value -1265.02
## 
## F = -1265.02
## final  value -1265.019277 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -385.81  |proj g|=           14
## At iterate    10  f =      -1070.3  |proj g|=        1.7684
## At iterate    20  f =      -1074.7  |proj g|=       0.64859
## ys=-1.518e-01  -gs= 1.371e-01, BFGS update SKIPPED
## ys=-7.603e+01  -gs= 3.615e+00, BFGS update SKIPPED
## At iterate    30  f =      -1177.2  |proj g|=        19.527
## At iterate    40  f =        -1262  |proj g|=      0.004444
## 
## iterations 41
## function evaluations 64
## segments explored during Cauchy searches 42
## BFGS updates skipped 2
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0044462
## final function value -1261.98
## 
## F = -1261.98
## final  value -1261.975372 
## converged

## 

## 
## ***************************

## Estimation of CV-set 9/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1305.7  |proj g|=       8.8586
## At iterate    10  f =        -1306  |proj g|=       0.26059
## 
## iterations 13
## function evaluations 26
## segments explored during Cauchy searches 13
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.051015
## final function value -1306.04
## 
## F = -1306.04
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1306.044405 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -396.77  |proj g|=           14
## At iterate    10  f =      -1105.6  |proj g|=        2.2523
## At iterate    20  f =      -1109.9  |proj g|=       0.25648
## ys=-6.297e-01  -gs= 1.942e+00, BFGS update SKIPPED
## At iterate    30  f =      -1302.6  |proj g|=       0.55208
## 
## iterations 35
## function evaluations 52
## segments explored during Cauchy searches 36
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00787603
## final function value -1302.6
## 
## F = -1302.6
## final  value -1302.604310 
## converged

## 

## 
## ***************************

## Estimation of CV-set 10/10

## Optimisation using starting value 1/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1292.3  |proj g|=        10.24
## At iterate    10  f =      -1293.8  |proj g|=       0.61724
## At iterate    20  f =      -1293.8  |proj g|=      0.020958
## 
## iterations 20
## function evaluations 41
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0209579
## final function value -1293.76
## 
## F = -1293.76
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1293.755426 
## converged

## 

## Optimisation using starting value 2/2

## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -387.69  |proj g|=           14
## At iterate    10  f =      -1081.4  |proj g|=        2.9417
## At iterate    20  f =      -1085.5  |proj g|=       0.22294
## 
## iterations 23
## function evaluations 37
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0124233
## final function value -1085.47
## 
## F = -1085.47
## final  value -1085.473201 
## converged

## 

print(est.denver.A.cv)

## Cross-validation parameter estimation for STmodel
##   with 10 CV-groups and 2 starting points.
##   Results: 10 converged, 0 not converged.
## 
## No fixed parameters.
## 
## Estimated function values and convergence info:
##       value convergence conv eigen.min eigen.all.min
## 1  1279.917        TRUE TRUE 1.8351328            NA
## 2  1280.933        TRUE TRUE 2.7925040            NA
## 3  1285.939        TRUE TRUE 2.9117065            NA
## 4  1274.043        TRUE TRUE 1.4687175            NA
## 5  1301.239        TRUE TRUE 1.3544770            NA
## 6  1273.617        TRUE TRUE 1.7810323            NA
## 7  1282.788        TRUE TRUE 1.3687804            NA
## 8  1265.019        TRUE TRUE 0.2122064            NA
## 9  1306.044        TRUE TRUE 1.8032551            NA
## 10 1293.755        TRUE TRUE 3.1516006            NA

par(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.A, pars="all"),
     plotCI((1:length(par))+.3, par, uiw=1.96*sd,
             col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.A.cv, "all", boxwex=.4, col="grey", add=TRUE)

#' Save the results as an .rdata object
save(denver.data.A, denver.model.A, est.denver.model.A, est.denver.A.cv,
     file = here::here("Results", "Denver_ST_Model_A.rdata"))

3.1.4 Prediction using the CV model

Making predictions using the CV model. Printing out the CV summary statistics as well

pred.A.cv <- predictCV(denver.model.A, est.denver.A.cv, LTA = T)
pred.A.cv.log <- predictCV(denver.model.A, est.denver.A.cv, LTA = T,
                           transform="unbiased")

head(pred.A.cv$pred.all$EX)

##            central         d_1        d_10         d_11        d_12        d_13
## 2008-12-29      NA -0.42509020 -0.29919410 -0.253778802 -0.29997252 -0.28480778
## 2009-01-05      NA -0.33745401 -0.22355879 -0.176026184 -0.22340549 -0.20997248
## 2009-01-12      NA -0.16860340 -0.06722598 -0.017699541 -0.07023096 -0.05849849
## 2009-01-19      NA -0.19323774 -0.10221236 -0.050581637 -0.09906798 -0.08896623
## 2009-01-26      NA -0.12695671 -0.04638523  0.007117754 -0.04189279 -0.03331327
## 2009-02-02      NA -0.03737667  0.03344767  0.088619327  0.03771543  0.04491357
##                   d_14         d_15      d_16         d_17       d_18
## 2008-12-29 -0.46606666 -0.388658381 0.2034593 -0.398988727 -0.4093148
## 2009-01-05 -0.38294825 -0.298634900 0.2514619 -0.309999522 -0.3254962
## 2009-01-12 -0.22334034 -0.124053289 0.3765865 -0.141111593 -0.1605016
## 2009-01-19 -0.24600816 -0.150822997 0.3208514 -0.162800490 -0.1885920
## 2009-01-26 -0.18307413 -0.082883369 0.3529177 -0.095162411 -0.1256528
## 2009-02-02 -0.09823997  0.009405756 0.4097405 -0.004809943 -0.0391435
##                   d_19       d_2         d_20        d_21        d_22
## 2008-12-29 -0.26739332 0.2419826 -0.233183779 -0.38157877 -0.43694213
## 2009-01-05 -0.18619766 0.2934278 -0.159571444 -0.29682928 -0.34550831
## 2009-01-12 -0.01872363 0.4280098 -0.005310516 -0.13154643 -0.17165784
## 2009-01-19 -0.05582689 0.3676511 -0.042091373 -0.15795028 -0.19436063
## 2009-01-26  0.00415006 0.4019541  0.011968596 -0.09411154 -0.12490642
## 2009-02-02  0.08976834 0.4629670  0.090167521 -0.00700837 -0.03199512
##                   d_23         d_25        d_26        d_28       d_29
## 2008-12-29 -0.03727008 -0.322228814 -0.31186712 -0.40101437 -0.5158686
## 2009-01-05  0.02486235 -0.246537633 -0.23289928 -0.32203488 -0.4293773
## 2009-01-12  0.16386225 -0.086030188 -0.07129028 -0.15829835 -0.2623840
## 2009-01-19  0.12143267 -0.126296122 -0.10573185 -0.19546780 -0.2871477
## 2009-01-26  0.16591946 -0.070954973 -0.04723558 -0.13723614 -0.2217779
## 2009-02-02  0.23401340  0.009901469  0.03573177 -0.05375665 -0.1332853
##                    d_3        d_31       d_32        d_33         d_34
## 2008-12-29 -0.30794846 -0.13373488 -0.5138657 -0.43435799 -0.236532203
## 2009-01-05 -0.22897077 -0.06982084 -0.4273476 -0.34848710 -0.162280634
## 2009-01-12 -0.07342900  0.07573486 -0.2597022 -0.17642206 -0.001625635
## 2009-01-19 -0.09999600  0.02876220 -0.2852506 -0.20912282 -0.045267901
## 2009-01-26 -0.04070177  0.07419080 -0.2199385 -0.14503622  0.008605008
## 2009-02-02  0.04082940  0.14484800 -0.1312758 -0.05568857  0.088684084
##                   d_35        d_36        d_37        d_38        d_39
## 2008-12-29 -0.40281945 -0.38209028 -0.38345643 -0.26410649 -0.32678654
## 2009-01-05 -0.31757678 -0.30158527 -0.30063421 -0.18952899 -0.25676739
## 2009-01-12 -0.15236790 -0.13984462 -0.13732926 -0.03360909 -0.10521652
## 2009-01-19 -0.17758478 -0.17105538 -0.16543749 -0.07040157 -0.14644021
## 2009-01-26 -0.11324000 -0.11102894 -0.10328177 -0.01558557 -0.09564504
## 2009-02-02 -0.02587609 -0.02716283 -0.01773624  0.06355225 -0.02011787
##                     d_4         d_40        d_41        d_42        d_43
## 2008-12-29 -0.280539857 -0.378186634 -0.36906123 -0.33153862 -0.34417146
## 2009-01-05 -0.193835910 -0.293015403 -0.29328637 -0.24724464 -0.26981802
## 2009-01-12  0.002584089 -0.098100517 -0.13690188 -0.05319120 -0.11411814
## 2009-01-19 -0.059795882 -0.161923794 -0.17164637 -0.11784055 -0.15112161
## 2009-01-26  0.001965910 -0.101509302 -0.11568549 -0.05819719 -0.09650256
## 2009-02-02  0.097571091 -0.007126749 -0.03576155  0.03548559 -0.01754347
##                   d_44         d_45        d_46        d_47         d_48
## 2008-12-29 -0.31504188 -0.306773132 -0.33081533 -0.34217411 -0.351014741
## 2009-01-05 -0.23591917 -0.236382577 -0.25049440 -0.25879469 -0.270463965
## 2009-01-12 -0.07415809 -0.085285401 -0.08544059 -0.08917624 -0.105184448
## 2009-01-19 -0.10845384 -0.125100094 -0.12134685 -0.12422313 -0.140874272
## 2009-01-26 -0.04982144 -0.073872157 -0.06193603 -0.06232659 -0.081261407
## 2009-02-02  0.03326945  0.001756799  0.02261350  0.02503365  0.003471469
##                    d_49         d_5        d_50        d_51        d_52
## 2008-12-29 -0.218660197 -0.34347179 -0.31382418 -0.36837408 -0.28965622
## 2009-01-05 -0.149111905 -0.26070363 -0.23805546 -0.28757662 -0.21187081
## 2009-01-12 -0.002829707 -0.09663337 -0.07747187 -0.12673288 -0.05351196
## 2009-01-19 -0.038276033 -0.12585175 -0.11766479 -0.15613565 -0.08636317
## 2009-01-26  0.012729499 -0.06384986 -0.06225548 -0.09569833 -0.02863494
## 2009-02-02  0.086738962  0.02184702  0.01866282 -0.01188030  0.05289279
##                    d_53        d_54        d_55        d_56        d_57
## 2008-12-29 -0.283092116 -0.32677770 -0.45949961 -0.33793421 -0.39809675
## 2009-01-05 -0.199533249 -0.24646550 -0.37284749 -0.25636081 -0.31974754
## 2009-01-12 -0.006201686 -0.07985904 -0.17647839 -0.09346380 -0.16074970
## 2009-01-19 -0.071543269 -0.11779423 -0.23890717 -0.12380726 -0.19318048
## 2009-01-26 -0.012546088 -0.05859387 -0.17719094 -0.06285560 -0.13496776
## 2009-02-02  0.080550305  0.02631969 -0.08162711  0.02188823 -0.05297001
##                   d_58          d_6        d_60         d_61        d_62
## 2008-12-29 -0.29602902 -0.341605339 -0.34345968 -0.371973834 -0.37267268
## 2009-01-05 -0.21791059 -0.262766807 -0.26210530 -0.286281982 -0.28595053
## 2009-01-12 -0.05913938 -0.101284785 -0.09447542 -0.120073744 -0.11672687
## 2009-01-19 -0.09178749 -0.135848126 -0.13142922 -0.145590216 -0.14386649
## 2009-01-26 -0.03377763 -0.077465522 -0.07131275 -0.080923111 -0.07855397
## 2009-02-02  0.04803602  0.005398679  0.01443214  0.006931764  0.01059889
##                   d_63       d_64        d_65        d_66         d_67
## 2008-12-29 -0.37681039 -0.4959411 -0.31895041 -0.33088536 -0.191873712
## 2009-01-05 -0.29518364 -0.4122248 -0.23832811 -0.25300314 -0.116906363
## 2009-01-12 -0.13352556 -0.2480418 -0.07647227 -0.09373081  0.067988421
## 2009-01-19 -0.16214743 -0.2753081 -0.10757258 -0.12755011 -0.005443494
## 2009-01-26 -0.10098115 -0.2123664 -0.04744304 -0.06984309  0.046001510
## 2009-02-02 -0.01650161 -0.1261076  0.03651665  0.01195635  0.132244565
##                   d_68          d_7         d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 -0.36359534 -0.380082266 -0.31851211   NA   NA   NA   NA  NA
## 2009-01-05 -0.27708875 -0.295911041 -0.24011260   NA   NA   NA   NA  NA
## 2009-01-12 -0.11063872 -0.131196041 -0.08162365   NA   NA   NA   NA  NA
## 2009-01-19 -0.13466541 -0.158144434 -0.11328450   NA   NA   NA   NA  NA
## 2009-01-26 -0.06920962 -0.094814001 -0.05495505   NA   NA   NA   NA  NA
## 2009-02-02  0.01916250 -0.008172109  0.02695015   NA   NA   NA   NA  NA

tail(pred.A.cv$pred.all$EX)

##            central       d_1      d_10      d_11      d_12      d_13      d_14
## 2020-11-02      NA 1.0507262 1.0086899 1.0671375 1.0404536 1.0471231 0.9267457
## 2020-11-09      NA 0.6495543 0.6084700 0.6480466 0.6275711 0.6160250 0.5264845
## 2020-11-16      NA 0.9930112 0.9705669 1.0050999 0.9899644 0.9873549 0.8734524
## 2020-11-23      NA 0.9106816 0.8809822 0.9241371 0.9018572 0.9019613 0.7762845
## 2020-11-30      NA 1.0528323 1.0259840 1.0789589 1.0540887 1.0631642 0.9233131
## 2020-12-07      NA        NA        NA        NA        NA        NA        NA
##                 d_15      d_16      d_17      d_18      d_19       d_2
## 2020-11-02 1.1210214 1.1371417 1.1236092 1.0229126 1.0957606 1.3054840
## 2020-11-09 0.7265924 0.7185737 0.7140727 0.6159464 0.7033632 0.8580272
## 2020-11-16 1.0895570 1.1127847 1.0676951 0.9712165 1.0664872 1.2371962
## 2020-11-23 0.9925457 1.0048637 0.9794156 0.8831923 0.9752908 1.1592842
## 2020-11-30 1.1249627 1.1944760 1.1244620 1.0297807 1.1122199 1.3150344
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_20      d_21      d_22      d_23      d_25      d_26
## 2020-11-02 1.0059040 1.0746504 1.1082391 1.0658122 0.9698567 1.0602988
## 2020-11-09 0.6064795 0.6652905 0.6985243 0.6689346 0.5706359 0.6294221
## 2020-11-16 0.9776844 1.0217154 1.0483121 1.0459459 0.9323433 0.9981913
## 2020-11-23 0.8909408 0.9347339 0.9629368 0.9524129 0.8412692 0.9165426
## 2020-11-30 1.0582063 1.0803400 1.1069639 1.0993921 0.9905781 1.0734375
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_28      d_29       d_3      d_31      d_32      d_33
## 2020-11-02 0.9248176 0.9534998 1.0844677 0.9818748 0.9623272 0.9882907
## 2020-11-09 0.5331957 0.5551532 0.6645467 0.5863394 0.5569678 0.5951057
## 2020-11-16 0.8833477 0.9043996 1.0318881 0.9542272 0.9117429 0.9368148
## 2020-11-23 0.7952065 0.8149507 0.9422594 0.8663463 0.8209447 0.8550988
## 2020-11-30 0.9398868 0.9579412 1.0953930 1.0129969 0.9667554 0.9939886
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_34      d_35      d_36      d_37      d_38      d_39
## 2020-11-02 1.0294080 1.0642775 0.9974888 1.0412300 1.0219983 0.8545244
## 2020-11-09 0.6253148 0.6541375 0.5928055 0.6345602 0.6173857 0.4726909
## 2020-11-16 0.9879691 1.0137144 0.9496459 0.9941513 0.9793151 0.8159157
## 2020-11-23 0.9019640 0.9213251 0.8613282 0.9032006 0.8882420 0.7353018
## 2020-11-30 1.0538837 1.0726021 1.0084492 1.0502736 1.0412946 0.8766090
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                  d_4      d_40      d_41      d_42      d_43      d_44
## 2020-11-02 1.0681120 0.7780578 0.9466048 0.9926037 0.9394276 1.0345445
## 2020-11-09 0.6580025 0.4083154 0.5407011 0.5885471 0.5380963 0.6262167
## 2020-11-16 1.0092313 0.7225507 0.9058082 0.9440133 0.8993429 0.9819527
## 2020-11-23 0.9395216 0.6279406 0.8108408 0.8694003 0.8076413 0.8994283
## 2020-11-30 1.0813790 0.7542407 0.9657512 1.0090471 0.9590759 1.0472659
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_45      d_46      d_47      d_48      d_49       d_5
## 2020-11-02 0.9096160 1.0072557 1.0539982 1.0031619 1.0166299 1.0012875
## 2020-11-09 0.5115848 0.6123676 0.6645520 0.6063096 0.6102352 0.6105413
## 2020-11-16 0.8716326 0.9556017 1.0225498 0.9557000 0.9820844 0.9428539
## 2020-11-23 0.7825943 0.8774532 0.9356096 0.8730891 0.8911619 0.8586572
## 2020-11-30 0.9329601 1.0189961 1.0618455 1.0152990 1.0380877 0.9989886
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_50      d_51      d_52      d_53      d_54      d_55
## 2020-11-02 0.9410239 1.0137597 1.0508521 1.1625364 1.0149022 0.8914866
## 2020-11-09 0.5581723 0.6196333 0.6393464 0.7178098 0.6139069 0.4989899
## 2020-11-16 0.8968702 0.9809427 0.9984167 1.0876010 0.9784398 0.8483559
## 2020-11-23 0.8192968 0.8897444 0.9143677 1.0129602 0.8936533 0.7700250
## 2020-11-30 0.9576246 1.0240487 1.0644770 1.1564030 1.0519860 0.9058540
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_56      d_57      d_58       d_6      d_60      d_61
## 2020-11-02 1.0128975 0.9317901 1.0767284 1.0112969 1.0041229 1.0596957
## 2020-11-09 0.6145393 0.5403825 0.6649401 0.6046561 0.6029944 0.6521355
## 2020-11-16 0.9535509 0.8878284 1.0268115 0.9653383 0.9527762 0.9958379
## 2020-11-23 0.8740066 0.8018455 0.9409341 0.8776155 0.8655090 0.9144330
## 2020-11-30 1.0179338 0.9448624 1.0858201 1.0249889 1.0113756 1.0610470
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_62      d_63      d_64      d_65      d_66      d_67
## 2020-11-02 1.0578885 1.0221486 0.9190848 1.0655690 1.0020878 1.0011043
## 2020-11-09 0.6563047 0.6177055 0.5248825 0.6279024 0.5975990 0.6023341
## 2020-11-16 0.9949120 0.9714508 0.8706186 0.9945009 0.9557804 0.9743314
## 2020-11-23 0.9145625 0.8853176 0.7833423 0.9168700 0.8659570 0.8928677
## 2020-11-30 1.0577934 1.0312216 0.9255889 1.0788609 1.0167814 1.0280331
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_68       d_7       d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.1669645 1.0304772 1.0327063   NA   NA   NA   NA  NA
## 2020-11-09 0.7444754 0.6205897 0.6260371   NA   NA   NA   NA  NA
## 2020-11-16 1.1049439 0.9682406 0.9829630   NA   NA   NA   NA  NA
## 2020-11-23 1.0173789 0.8835628 0.8977986   NA   NA   NA   NA  NA
## 2020-11-30 1.1628995 1.0319215 1.0454443   NA   NA   NA   NA  NA
## 2020-12-07        NA        NA        NA   NA   NA   NA   NA  NA

head(pred.A.cv.log$pred.all$EX)

##            central       d_1      d_10      d_11      d_12      d_13      d_14
## 2008-12-29      NA 0.6706146 0.7606878 0.7965575 0.7612359 0.7728678 0.6447413
## 2009-01-05      NA 0.7318498 0.8202008 0.8606169 0.8214333 0.8325421 0.7002994
## 2009-01-12      NA 0.8662736 0.9587303 1.0078974 0.9570165 0.9683108 0.8211545
## 2009-01-19      NA 0.8450322 0.9255531 0.9750002 0.9294927 0.9389298 0.8024734
## 2009-01-26      NA 0.9027973 0.9785070 1.0326554 0.9839063 0.9923840 0.8543571
## 2009-02-02      NA 0.9872792 1.0596663 1.1201223 1.0651955 1.0728906 0.9297896
##                 d_15     d_16      d_17      d_18      d_19      d_2      d_20
## 2008-12-29 0.6946862 1.259380 0.6883669 0.6813451 0.7848421 1.309593 0.8131327
## 2009-01-05 0.7600640 1.320701 0.7522142 0.7406528 0.8510093 1.378013 0.8748953
## 2009-01-12 0.9049800 1.496132 0.8903847 0.8732426 1.0059346 1.575808 1.0204620
## 2009-01-19 0.8810218 1.414538 0.8710917 0.8488276 0.9691096 1.482930 0.9833135
## 2009-01-26 0.9429111 1.460218 0.9318822 0.9037710 1.0288504 1.534191 1.0376768
## 2009-02-02 1.0340324 1.545246 1.0198560 0.9852654 1.1206796 1.630296 1.1218578
##                 d_21      d_22      d_23      d_25      d_26      d_28
## 2008-12-29 0.7005308 0.6641760 0.9899413 0.7424012 0.7526664 0.6861555
## 2009-01-05 0.7622531 0.7273895 1.0529158 0.8007109 0.8140902 0.7424850
## 2009-01-12 0.8990056 0.8651092 1.2094413 0.9400527 0.9564450 0.8745136
## 2009-01-19 0.8753760 0.8453610 1.1587992 0.9028977 0.9237048 0.8425539
## 2009-01-26 0.9329034 0.9058720 1.2111710 0.9542259 0.9790364 0.8930292
## 2009-02-02 1.0176520 0.9938189 1.2962252 1.0345452 1.0634588 0.9707403
##                 d_29       d_3      d_31      d_32      d_33      d_34
## 2008-12-29 0.6124998 0.7551884 0.8974433 0.6137073 0.6641565 0.8094408
## 2009-01-05 0.6676277 0.8168745 0.9564308 0.6689308 0.7235239 0.8716082
## 2009-01-12 0.7887518 0.9539608 1.1060374 0.7907745 0.8591761 1.0232820
## 2009-01-19 0.7692809 0.9286305 1.0550829 0.7706217 0.8313757 0.9793967
## 2009-01-26 0.8210921 0.9850788 1.1039452 0.8224522 0.8862606 1.0334441
## 2009-02-02 0.8969287 1.0685177 1.1846203 0.8985467 0.9689700 1.1194652
##                 d_35      d_36      d_37      d_38      d_39       d_4
## 2008-12-29 0.6857350 0.7001492 0.6996790 0.7878260 0.7398876 0.7752066
## 2009-01-05 0.7465360 0.7585759 0.7597895 0.8485247 0.7933479 0.8451613
## 2009-01-12 0.8804185 0.8914687 0.8942559 0.9913879 0.9229616 1.0283162
## 2009-01-19 0.8583079 0.8638445 0.8692075 0.9553203 0.8855175 0.9659078
## 2009-01-26 0.9151873 0.9170848 0.9247184 1.0089276 0.9315126 1.0272508
## 2009-02-02 0.9985963 0.9971405 1.0071074 1.0918260 1.0044663 1.1301398
##                 d_40      d_41      d_42      d_43      d_44      d_45
## 2008-12-29 0.7030886 0.7098239 0.7366633 0.7272078 0.7502807 0.7554435
## 2009-01-05 0.7653613 0.7653929 0.8012066 0.7830606 0.8116354 0.8102097
## 2009-01-12 0.9298222 0.8946382 0.9725317 0.9147006 0.9537060 0.9420288
## 2009-01-19 0.8721318 0.8638274 0.9114382 0.8812370 0.9211939 0.9049858
## 2009-01-26 0.9262704 0.9133193 0.9672706 0.9305039 0.9765080 0.9523179
## 2009-02-02 1.0178001 0.9891166 1.0621082 1.0067787 1.0608435 1.0269316
##                 d_46      d_47      d_48      d_49       d_5      d_50
## 2008-12-29 0.7360538 0.7282917 0.7213351 0.8257199 0.7276448 0.7486671
## 2009-01-05 0.7975489 0.7914178 0.7817802 0.8847847 0.7902313 0.8075315
## 2009-01-12 0.9406072 0.9375028 0.9222181 1.0237441 0.9309176 0.9481325
## 2009-01-19 0.9073775 0.9050421 0.8898306 0.9877512 0.9039379 0.9107246
## 2009-01-26 0.9628709 0.9626796 0.9444417 1.0391444 0.9616061 0.9625634
## 2009-02-02 1.0477803 1.0504310 1.0279144 1.1187164 1.0475156 1.0436490
##                 d_51      d_52      d_53      d_54      d_55      d_56
## 2008-12-29 0.7097669 0.7684857 0.7732306 0.7395916 0.6481809 0.7316854
## 2009-01-05 0.7692716 0.8303148 0.8403598 0.8012357 0.7066362 0.7936705
## 2009-01-12 0.9032798 0.9724409 1.0193212 0.9462785 0.8597275 0.9338729
## 2009-01-19 0.8769167 0.9407300 0.9546273 0.9108793 0.8075114 0.9057879
## 2009-01-26 0.9313829 0.9963874 1.0124510 0.9662798 0.8587558 0.9625627
## 2009-02-02 1.0126707 1.0808106 1.1110667 1.0517828 0.9447293 1.0475588
##                 d_57      d_58       d_6      d_60      d_61      d_62
## 2008-12-29 0.6890545 0.7630993 0.7306130 0.7273561 0.7072918 0.7082638
## 2009-01-05 0.7449823 0.8248466 0.7901349 0.7888020 0.7703354 0.7720273
## 2009-01-12 0.8731314 0.9665146 0.9281829 0.9325480 0.9093790 0.9139599
## 2009-01-19 0.8450733 0.9352523 0.8963010 0.8985438 0.8862629 0.8891429
## 2009-01-26 0.8955566 0.9909217 0.9498831 0.9540676 0.9452884 0.9488497
## 2009-02-02 0.9719376 1.0752384 1.0316852 1.0393537 1.0319375 1.0370641
##                 d_63      d_64      d_65      d_66      d_67      d_68
## 2008-12-29 0.7038043 0.6252409 0.7457820 0.7368612 0.8470805 0.7131668
## 2009-01-05 0.7634419 0.6795636 0.8081115 0.7963399 0.9127455 0.7773820
## 2009-01-12 0.8971649 0.8005343 0.9497917 0.9336236 1.0978207 0.9179349
## 2009-01-19 0.8716607 0.7787665 0.9204619 0.9024040 1.0198602 0.8959479
## 2009-01-26 0.9264755 0.8291530 0.9772923 0.9558602 1.0734972 0.9563837
## 2009-02-02 1.0080016 0.9036717 1.0627033 1.0372061 1.1700129 1.0445998
##                  d_7       d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.7015800 0.7460545   NA   NA   NA   NA  NA
## 2009-01-05 0.7629534 0.8066647   NA   NA   NA   NA  NA
## 2009-01-12 0.8993207 0.9449591   NA   NA   NA   NA  NA
## 2009-01-19 0.8752061 0.9153104   NA   NA   NA   NA  NA
## 2009-01-26 0.9322483 0.9701141   NA   NA   NA   NA  NA
## 2009-02-02 1.0164684 1.0527666   NA   NA   NA   NA  NA

tail(pred.A.cv.log$pred.all$EX)

##            central      d_1     d_10     d_11     d_12     d_13     d_14
## 2020-11-02      NA 2.882844 2.758274 2.927953 2.849952 2.866982 2.543824
## 2020-11-09      NA 1.930042 1.848378 1.925364 1.885719 1.862744 1.704543
## 2020-11-16      NA 2.720860 2.654734 2.751349 2.709094 2.700110 2.411326
## 2020-11-23      NA 2.505726 2.427134 2.537214 2.480444 2.478935 2.187893
## 2020-11-30      NA 2.888383 2.805759 2.961924 2.888140 2.912394 2.534285
## 2020-12-07      NA       NA       NA       NA       NA       NA       NA
##                d_15     d_16     d_17     d_18     d_19      d_2     d_20
## 2020-11-02 3.082275 3.132966 3.096215 2.798747 3.006990 3.711721 2.749532
## 2020-11-09 2.077604 2.061225 2.055608 1.862871 2.030897 2.372429 1.843946
## 2020-11-16 2.986676 3.056966 2.927466 2.657325 2.919903 3.465943 2.672555
## 2020-11-23 2.710504 2.744040 2.679987 2.433282 2.665290 3.205905 2.450346
## 2020-11-30 3.094232 3.316756 3.098202 2.817318 3.056313 3.745974 2.896332
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_21     d_22     d_23     d_25     d_26     d_28     d_29
## 2020-11-02 2.948014 3.055046 2.918045 2.652956 2.907973 2.538925 2.613302
## 2020-11-09 1.957553 2.027810 1.961922 1.779674 1.889771 1.716171 1.754513
## 2020-11-16 2.795636 2.876708 2.860074 2.555170 2.732254 2.435688 2.487747
## 2020-11-23 2.562619 2.641095 2.604509 2.332708 2.517836 2.230161 2.274773
## 2020-11-30 2.964170 3.050051 3.016710 2.708316 2.945363 2.577299 2.624347
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_3     d_31     d_32     d_33     d_34     d_35     d_36
## 2020-11-02 2.976240 2.684562 2.633946 2.708008 2.815652 2.915474 2.728219
## 2020-11-09 1.955466 1.807451 1.755997 1.827529 1.879557 1.934444 1.820077
## 2020-11-16 2.823231 2.611054 2.503633 2.571844 2.701046 2.771363 2.600361
## 2020-11-23 2.581016 2.391288 2.286197 2.369946 2.478351 2.526675 2.380422
## 2020-11-30 3.007955 2.768900 2.644962 2.722974 2.884871 2.939221 2.757583
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_37     d_38     d_39      d_4     d_40     d_41     d_42
## 2020-11-02 2.850910 2.794233 2.370156 2.934581 2.193477 2.592523 2.719307
## 2020-11-09 1.898133 1.864248 1.617786 1.947176 1.515386 1.727421 1.815290
## 2020-11-16 2.719326 2.677057 2.280132 2.766407 2.074761 2.488447 2.589984
## 2020-11-23 2.482762 2.443888 2.103451 2.580007 1.887377 2.262860 2.403657
## 2020-11-30 2.875985 2.847947 2.422631 2.973119 2.141378 2.641879 2.763781
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_43     d_44     d_45     d_46     d_47     d_48     d_49
## 2020-11-02 2.573452 2.837354 2.500204 2.760023 2.884754 2.746439 2.779292
## 2020-11-09 1.722592 1.885929 1.679076 1.859538 1.954098 1.846754 1.850934
## 2020-11-16 2.471951 2.691393 2.406597 2.620969 2.795120 2.619024 2.684384
## 2020-11-23 2.255228 2.478012 2.201445 2.423908 2.562271 2.411325 2.450908
## 2020-11-30 2.623847 2.872641 2.558523 2.792431 2.906929 2.779788 2.838648
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_5     d_50     d_51     d_52     d_53     d_54     d_55
## 2020-11-02 2.744472 2.584436 2.770826 2.879907 3.216428 2.774320 2.458789
## 2020-11-09 1.856660 1.762323 1.868146 1.908189 2.061572 1.857712 1.660464
## 2020-11-16 2.588403 2.472705 2.681023 2.732310 2.983805 2.674674 2.354676
## 2020-11-23 2.379295 2.288108 2.447224 2.511893 2.769067 2.457145 2.177168
## 2020-11-30 2.737660 2.627523 2.798887 2.918579 3.196047 2.878588 2.493821
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_56     d_57     d_58      d_6     d_60     d_61     d_62
## 2020-11-02 2.776892 2.558377 2.952317 2.769914 2.747011 2.908258 2.908396
## 2020-11-09 1.864346 1.729599 1.955655 1.844211 1.839182 1.934634 1.946232
## 2020-11-16 2.616589 2.448009 2.808177 2.644907 2.609221 2.727976 2.730283
## 2020-11-23 2.416420 2.246209 2.576958 2.422584 2.391077 2.514584 2.519293
## 2020-11-30 2.790393 2.591462 2.978609 2.807082 2.766486 2.911542 2.907074
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_63     d_64     d_65     d_66     d_67     d_68      d_7
## 2020-11-02 2.798118 2.526940 2.920288 2.740792 2.736374 3.230491 2.822425
## 2020-11-09 1.867185 1.703542 1.885002 1.828861 1.836364 2.117154 1.873171
## 2020-11-16 2.659452 2.406970 2.719528 2.616464 2.663721 3.035826 2.651759
## 2020-11-23 2.439862 2.205669 2.516257 2.391593 2.455210 2.781173 2.436340
## 2020-11-30 2.823025 2.542705 2.958610 2.780836 2.810434 3.216702 2.825872
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.827221   NA   NA   NA   NA  NA
## 2020-11-09 1.882410   NA   NA   NA   NA  NA
## 2020-11-16 2.689678   NA   NA   NA   NA  NA
## 2020-11-23 2.469984   NA   NA   NA   NA  NA
## 2020-11-30 2.862860   NA   NA   NA   NA  NA
## 2020-12-07       NA   NA   NA   NA   NA  NA

summary(pred.A.cv)

## Cross-validation predictions for STmodel with 10 CV-groups.
##   Predictions for 712 observations.
##   Temporal averages for 63 locations.
## 
## RMSE:
##              EX.mu EX.mu.beta         EX
## obs     0.16720556 0.16720556 0.11062431
## average 0.09258214 0.09258214 0.04833502
## 
## R2:
##             EX.mu EX.mu.beta        EX
## obs     0.5990090  0.5990090 0.8244769
## average 0.6215252  0.6215252 0.8968412
## 
## Coverage of 95% prediction intervals:
##                EX
## obs     0.9536517
## average 0.9047619

summary(pred.A.cv.log)

## Cross-validation predictions for STmodel with 10 CV-groups.
##   Predictions for 712 observations.
##   Temporal averages for 63 locations.
## 
## RMSE:
##             EX.mu EX.mu.beta         EX    EX.pred
## obs     0.2211963  0.2211963 0.15552805 0.15542748
## average 0.1045434  0.1045434 0.05577006 0.05597244
## 
## R2:
##             EX.mu EX.mu.beta        EX   EX.pred
## obs     0.6185378  0.6185378 0.8114124 0.8116562
## average 0.6822837  0.6822837 0.9095833 0.9089259
## 
## Coverage of 95% prediction intervals:
##           EX.pred
## obs     0.9536517
## average 0.9523810

par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
     col=c("ID", "black", "grey"),
     ylim=c(-1,2),
     xlab="Observations", ylab="Predictions",
     main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
                                    xlab="Observations", ylab="Predictions",
                                    main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")

jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModA.jpeg"),
     width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
     col=c("ID", "black", "grey"),
     ylim=c(-1,2),
     xlab="Observations", ylab="Predictions",
     main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
                                      xlab="Observations", ylab="Predictions",
                                      main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()

## quartz_off_screen 
##                 2

3.1.5 Predictions for 2009-2020

What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites

x.A <- coef(est.denver.model.A, pars = "cov")$par
x.A

## [1] -7.558680 -7.132387 12.520027 -3.032854 -4.759856

E.A <- predict(denver.model.A, est.denver.model.A, LTA = T,
               transform="unbiased", pred.var=FALSE)
print(E.A)

## Prediction for STmodel.
## 
## Regression parameters:
##   0 Spatio-temporal covariate(s).
##   16 beta-fields regression parameters in x$pars.
## 
## Regression parameters are assumed to be known and
##  prediction variances do NOT include
##  uncertainties in regression parameters.
## 
## Prediction of beta-fields, (x$beta):
## List of 2
##  $ mu: num [1:69, 1:2] 0.2211736 0.0025188 -0.0001717 0.0546205 0.0000313 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX: num [1:69, 1:2] 0.23143 0.01497 -0.00252 0.03873 -0.00599 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Predictions for log-Gaussian field of type: unbiased 
## 
## Predictions for 624 times at 69 locations.
## List of 5
##  $ EX.mu     : num [1:624, 1:69] 1.19 1.25 1.43 1.35 1.39 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX.mu.beta: num [1:624, 1:69] 1.25 1.31 1.49 1.4 1.44 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX        : num [1:624, 1:69] 1.28 1.34 1.53 1.43 1.48 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX.pred   : num [1:624, 1:69] 1.29 1.35 1.53 1.44 1.48 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ log.EX    : num [1:624, 1:69] 0.223 0.27 0.399 0.336 0.367 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Variances have been computed.
## List of 2
##  $ log.VX     : num [1:624, 1:69] 0.0487 0.0486 0.0486 0.0485 0.0485 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ log.VX.pred: num [1:624, 1:69] 0.0573 0.0572 0.0571 0.0571 0.057 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Mean squared prediciton errors NOT been computed.
## 
## 69 temporal averages have been compute.
## List of 1
##  $ LTA:'data.frame': 69 obs. of  4 variables:

week_preds <- as.data.frame(E.A$EX) %>%
  mutate(week = as.Date(rownames(E.A$EX)),
         year = year(as.Date(rownames(E.A$EX))))

week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]

box_preds <- week_preds %>%
  select(week, all_of(week_sites)) %>%
  #filter(week %in% week_weeks) %>%
  mutate(month = month(week),
         year = year(week)) %>%
  filter(year %in% c(2009:2020))

box_data <- box_preds %>%
  pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)

##       week                month             year        location        
##  Min.   :2009-01-05   Min.   : 1.000   Min.   :2009   Length:42364      
##  1st Qu.:2011-12-26   1st Qu.: 4.000   1st Qu.:2011   Class :character  
##  Median :2014-12-22   Median : 7.000   Median :2014   Mode  :character  
##  Mean   :2014-12-22   Mean   : 6.506   Mean   :2014                     
##  3rd Qu.:2017-12-18   3rd Qu.: 9.000   3rd Qu.:2017                     
##  Max.   :2020-12-07   Max.   :12.000   Max.   :2020                     
##                                                                         
##       pred       
##  Min.   :0.3344  
##  1st Qu.:1.1149  
##  Median :1.3724  
##  Mean   :1.4213  
##  3rd Qu.:1.6580  
##  Max.   :3.8895  
##  NA's   :1080

#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
  group_by(month) %>%
  summarize(median = median(pred, na.rm=T)) %>%
  arrange(desc(median))
box_summary

pred_box_plot <- ggplot(box_data) +
  geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
               show.legend = F) +
  #scale_color_viridis(name = "Month", discrete = T) +
  facet_wrap(. ~ year, scales = "free") +
  xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
  # theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot

3.2 Model B: 2 temporal trend (basis) functions

3.2.1 Create the model object

For this version of the model, use iid for cov.beta (beta, beta1, beta2) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.

denver.data.B <- denver.data2.2
names(denver.data.B$covars)

##  [1] "ID"              "lon"             "lat"             "impervious_2500"
##  [5] "open_2500"       "med_int_2500"    "high_int_100"    "dist_m_cafos"   
##  [9] "aadt_100"        "aadt_250"        "x"               "y"              
## [13] "type"

LUR.B <- list(covar_fun, covar_fun, covar_fun)

cov.beta.B <-  list(covf = c("iid", "iid", "iid"), nugget = T)
cov.nu.B <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.B <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")

denver.model.B <- createSTmodel(denver.data.B, LUR = LUR.B,
                                ST = c("bc_st_no2", "bc_st_smk"),
                                cov.beta = cov.beta.B, cov.nu = cov.nu.B,
                                locations = locations.B)

## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.

denver.model.B

## STmodel-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 624 (observed: 231)
##  No. obs: 943
## 
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## Models for the beta-fields are:
## $const
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 + 
##     dist_m_cafos + aadt_100 + aadt_250
## 
## $V1
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 + 
##     dist_m_cafos + aadt_100 + aadt_250
## 
## $V2
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 + 
##     dist_m_cafos + aadt_100 + aadt_250
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## Covariance model for the beta-field(s):
##  Covariance type(s): iid, iid, iid 
##  Nugget: Yes, Yes, Yes 
## Covariance model for the nu-field(s):
##  Covariance type: exp 
##  Nugget: ~1 
##  Random effect: No 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 712
##   Dates: 2018-05-07 to 2020-04-06

3.2.2 Estimate model parameters

• nugget values can range from -1 to -6
• range values can range from 0 to 4
• try starting values where the sill and nugget are the same (e.g., both 0) and where they are very different (e.g., 0 and -5)
• make sure the starting values are pretty different

names <- loglikeSTnames(denver.model.B, all=FALSE)
names

## [1] "log.nugget.const.iid"          "log.nugget.V1.iid"            
## [3] "log.nugget.V2.iid"             "nu.log.range.exp"             
## [5] "nu.log.sill.exp"               "nu.log.nugget.(Intercept).exp"

x.init.B <- cbind(c(0, 0, 0, 0, 0, 0),
                  c(-10, -5, -5, 8, -3, -5),
                  c(-10, -5, -5, 8, -5, -5),
                  c(-10, -5, -5, 10, -3, -5),
                  c(-10, -5, -5, 10, -5, -5),
                  c(-12, -5, -5, 8, -5, -5),
                  c(-12, -5, -5, 10, -5, -5),
                  c(-12, -5, -5, 12, -5, -5))

rownames(x.init.B) <- loglikeSTnames(denver.model.B, all=FALSE)
x.init.B

##                               [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## log.nugget.const.iid             0  -10  -10  -10  -10  -12  -12  -12
## log.nugget.V1.iid                0   -5   -5   -5   -5   -5   -5   -5
## log.nugget.V2.iid                0   -5   -5   -5   -5   -5   -5   -5
## nu.log.range.exp                 0    8    8   10   10    8   10   12
## nu.log.sill.exp                  0   -3   -5   -3   -5   -5   -5   -5
## nu.log.nugget.(Intercept).exp    0   -5   -5   -5   -5   -5   -5   -5

#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.B <- estimate.STmodel(denver.model.B, x.init.B)

## Optimisation using starting value 1/8

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=       474.52  |proj g|=           15
## At iterate    10  f =      -1257.9  |proj g|=       0.62776
## At iterate    20  f =      -1258.2  |proj g|=       0.34897
## 
## iterations 29
## function evaluations 35
## segments explored during Cauchy searches 34
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 2
## norm of the final projected gradient 0.000608333
## final function value -1258.28
## 
## F = -1258.28
## final  value -1258.282191 
## converged

## Optimisation using starting value 2/8

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1145.5  |proj g|=           12
## At iterate    10  f =      -1458.9  |proj g|=        5.2642
## At iterate    20  f =      -1459.4  |proj g|=      0.069819
## At iterate    30  f =      -1459.5  |proj g|=       0.85468
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1459.490471 
## stopped after 39 iterations

## Optimisation using starting value 3/8

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1143.1  |proj g|=           20
## At iterate    10  f =      -1456.1  |proj g|=         3.879
## At iterate    20  f =      -1459.4  |proj g|=       0.14248
## At iterate    30  f =      -1459.4  |proj g|=        1.1349
## At iterate    40  f =      -1459.5  |proj g|=      0.072186
## 
## iterations 48
## function evaluations 61
## segments explored during Cauchy searches 53
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0304898
## final function value -1459.49
## 
## F = -1459.49
## final  value -1459.490601 
## converged

## Optimisation using starting value 4/8

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1309.4  |proj g|=           12
## At iterate    10  f =      -1458.9  |proj g|=        11.518
## At iterate    20  f =      -1459.4  |proj g|=       0.96479
## At iterate    30  f =      -1459.5  |proj g|=      0.031888
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1459.490603 
## stopped after 30 iterations

## Optimisation using starting value 5/8

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1200.7  |proj g|=           20
## At iterate    10  f =      -1459.2  |proj g|=        1.2752
## At iterate    20  f =      -1459.4  |proj g|=      0.078369
## At iterate    30  f =      -1459.4  |proj g|=       0.56669
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1459.490611 
## stopped after 38 iterations

## Optimisation using starting value 6/8

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1142.8  |proj g|=           20
## At iterate    10  f =      -1456.6  |proj g|=        10.275
## At iterate    20  f =      -1459.3  |proj g|=      0.037973
## 
## iterations 25
## function evaluations 40
## segments explored during Cauchy searches 30
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0493346
## final function value -1459.35
## 
## F = -1459.35
## final  value -1459.348069 
## converged

## Optimisation using starting value 7/8

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1200.4  |proj g|=           20
## At iterate    10  f =      -1459.2  |proj g|=        1.2461
## At iterate    20  f =      -1459.3  |proj g|=      0.028556
## 
## iterations 20
## function evaluations 36
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0285564
## final function value -1459.35
## 
## F = -1459.35
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1459.347897 
## converged

## Optimisation using starting value 8/8

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=        -1200  |proj g|=           20
## At iterate    10  f =        -1459  |proj g|=         4.961
## At iterate    20  f =      -1459.3  |proj g|=      0.026597
## At iterate    30  f =      -1459.3  |proj g|=        0.4351
## ys=-7.527e-03  -gs= 1.809e-02, BFGS update SKIPPED
## At iterate    40  f =      -1459.5  |proj g|=       0.24613
## 
## iterations 46
## function evaluations 63
## segments explored during Cauchy searches 50
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0258296
## final function value -1459.49
## 
## F = -1459.49
## final  value -1459.490552 
## converged

print(est.denver.model.B)

## Optimisation for STmodel with 8 starting points.
##   Results: 2 converged, 6 not converged, 0 failed.
##   Best result for starting point 5, optimisation has NOT converged
##   Best converged result for starting point 3 
## 
## No fixed parameters.
## 
## Estimated parameters for all starting point(s):
##                                        [,1]          [,2]          [,3]
## gamma.bc_st_no2                 0.006825267  0.0086902950  0.0086810932
## gamma.bc_st_smk                 0.050977166  0.0490210509  0.0490347152
## alpha.const.(Intercept)         0.175407462  0.1956111940  0.1958073526
## alpha.const.impervious_2500     0.002096213  0.0152705785  0.0153226189
## alpha.const.open_2500           0.007918598  0.0064206957  0.0064102731
## alpha.const.med_int_2500        0.016339133 -0.0003250761 -0.0003660604
## alpha.const.high_int_100        0.002482419  0.0048685312  0.0048472442
## alpha.const.dist_m_cafos       -0.021717818 -0.0214981651 -0.0215274974
## alpha.const.aadt_100            0.030277467  0.0166623609  0.0166407242
## alpha.const.aadt_250            0.022896805  0.0178550584  0.0178817712
## alpha.V1.(Intercept)            0.166162945  0.1783673742  0.1783837079
## alpha.V1.impervious_2500        0.014938297  0.0075446851  0.0075584878
## alpha.V1.open_2500             -0.006466654  0.0006784633  0.0006804785
## alpha.V1.med_int_2500          -0.007280071 -0.0036866182 -0.0036926297
## alpha.V1.high_int_100           0.023599868  0.0238944621  0.0238902363
## alpha.V1.dist_m_cafos           0.004421040 -0.0024551014 -0.0024314856
## alpha.V1.aadt_100              -0.015712578 -0.0206936558 -0.0206953131
## alpha.V1.aadt_250              -0.031984221 -0.0200627761 -0.0200603014
## alpha.V2.(Intercept)            0.157566420  0.1723380449  0.1723940701
## alpha.V2.impervious_2500       -0.022409303 -0.0170009993 -0.0169686559
## alpha.V2.open_2500              0.013478857  0.0117515720  0.0117436817
## alpha.V2.med_int_2500           0.026519369  0.0178107329  0.0177856997
## alpha.V2.high_int_100          -0.006151440 -0.0041084710 -0.0041163634
## alpha.V2.dist_m_cafos           0.010196517  0.0100957490  0.0100791610
## alpha.V2.aadt_100              -0.021198667 -0.0268102887 -0.0268312866
## alpha.V2.aadt_250               0.004345172  0.0021175995  0.0021302273
## log.nugget.const.iid          -14.319024162 -8.7133077165 -8.6933320465
## log.nugget.V1.iid             -15.000000000 -7.0772541944 -7.0746701234
## log.nugget.V2.iid             -15.000000000 -7.5629631080 -7.5683449287
## nu.log.range.exp                0.000000000 13.8028190119 13.8120753074
## nu.log.sill.exp                -4.623411406 -3.0088155216 -3.0085177177
## nu.log.nugget.(Intercept).exp  -4.149419586 -4.8187117583 -4.8182150293
##                                        [,4]          [,5]           [,6]
## gamma.bc_st_no2                0.0086853418  0.0086804968   0.0089507205
## gamma.bc_st_smk                0.0490259149  0.0490321152   0.0490594641
## alpha.const.(Intercept)        0.1956978301  0.1958058811   0.1922766174
## alpha.const.impervious_2500    0.0153057073  0.0153252484   0.0140527012
## alpha.const.open_2500          0.0064092579  0.0064068806   0.0069662685
## alpha.const.med_int_2500      -0.0003602946 -0.0003741714   0.0012286750
## alpha.const.high_int_100       0.0048541069  0.0048462292   0.0053625754
## alpha.const.dist_m_cafos      -0.0215188945 -0.0215273574  -0.0209388328
## alpha.const.aadt_100           0.0166421541  0.0166374468   0.0174753357
## alpha.const.aadt_250           0.0178793006  0.0178854015   0.0169291157
## alpha.V1.(Intercept)           0.1783822229  0.1783840296   0.1777210937
## alpha.V1.impervious_2500       0.0075506172  0.0075523751   0.0077571785
## alpha.V1.open_2500             0.0006825803  0.0006810681   0.0005800503
## alpha.V1.med_int_2500         -0.0036849352 -0.0036871961  -0.0040671677
## alpha.V1.high_int_100          0.0238863602  0.0238876258   0.0242013537
## alpha.V1.dist_m_cafos         -0.0024432810 -0.0024368070  -0.0024536176
## alpha.V1.aadt_100             -0.0206952794 -0.0206947608  -0.0206231665
## alpha.V1.aadt_250             -0.0200593952 -0.0200573132  -0.0203550796
## alpha.V2.(Intercept)           0.1723653835  0.1723915816   0.1717389221
## alpha.V2.impervious_2500      -0.0169806291 -0.0169680165  -0.0176446886
## alpha.V2.open_2500             0.0117443581  0.0117413398   0.0121515502
## alpha.V2.med_int_2500          0.0177917264  0.0177816580   0.0186519515
## alpha.V2.high_int_100         -0.0041138614 -0.0041172903  -0.0038791942
## alpha.V2.dist_m_cafos          0.0100859408  0.0100788070   0.0104441570
## alpha.V2.aadt_100             -0.0268276554 -0.0268338080  -0.0261199038
## alpha.V2.aadt_250              0.0021296855  0.0021331199   0.0015899649
## log.nugget.const.iid          -8.6940671217 -8.6870154341 -11.9673583758
## log.nugget.V1.iid             -7.0766936210 -7.0748494524  -7.0723990876
## log.nugget.V2.iid             -7.5678727748 -7.5681824916  -7.4631928925
## nu.log.range.exp              13.8060928682 13.8095411299  13.8321924769
## nu.log.sill.exp               -3.0090198820 -3.0084829692  -3.0118696588
## nu.log.nugget.(Intercept).exp -4.8184536221 -4.8183782834  -4.8080020534
##                                         [,7]          [,8]
## gamma.bc_st_no2                 0.0089512000  0.0086805587
## gamma.bc_st_smk                 0.0490578896  0.0490258663
## alpha.const.(Intercept)         0.1922668708  0.1957613345
## alpha.const.impervious_2500     0.0140465974  0.0153274410
## alpha.const.open_2500           0.0069623102  0.0063987275
## alpha.const.med_int_2500        0.0012280383 -0.0003880943
## alpha.const.high_int_100        0.0053654093  0.0048444566
## alpha.const.dist_m_cafos       -0.0209381425 -0.0215300089
## alpha.const.aadt_100            0.0174719150  0.0166265606
## alpha.const.aadt_250            0.0169329549  0.0179017693
## alpha.V1.(Intercept)            0.1777247001  0.1784124588
## alpha.V1.impervious_2500        0.0077603418  0.0075465152
## alpha.V1.open_2500              0.0005855862  0.0006816656
## alpha.V1.med_int_2500          -0.0040655491 -0.0036764101
## alpha.V1.high_int_100           0.0241961668  0.0238800331
## alpha.V1.dist_m_cafos          -0.0024543853 -0.0024389340
## alpha.V1.aadt_100              -0.0206214482 -0.0207037662
## alpha.V1.aadt_250              -0.0203575991 -0.0200540593
## alpha.V2.(Intercept)            0.1717365101  0.1723802107
## alpha.V2.impervious_2500       -0.0176526009 -0.0169684392
## alpha.V2.open_2500              0.0121509475  0.0117355998
## alpha.V2.med_int_2500           0.0186569176  0.0177763912
## alpha.V2.high_int_100          -0.0038770257 -0.0041188040
## alpha.V2.dist_m_cafos           0.0104490010  0.0100798319
## alpha.V2.aadt_100              -0.0261205759 -0.0268426034
## alpha.V2.aadt_250               0.0015905254  0.0021414966
## log.nugget.const.iid          -11.9786223678 -8.6744466377
## log.nugget.V1.iid              -7.0738093944 -7.0784616664
## log.nugget.V2.iid              -7.4661100356 -7.5700023965
## nu.log.range.exp               13.8307608958 13.8069718049
## nu.log.sill.exp                -3.0117132930 -3.0090709173
## nu.log.nugget.(Intercept).exp  -4.8080881154 -4.8185034029
## 
## Function value(s):
## [1] 1258.282 1459.490 1459.491 1459.491 1459.491 1459.348 1459.348 1459.491

3.2.3 Cross-validation

Define the CV groups

set.seed(1000)
site_idsB <- unique(denver.model.B$obs$ID[which(str_detect(denver.model.B$obs$ID, "d_") == T)])
site_idsB

##  [1] "d_1"  "d_18" "d_22" "d_25" "d_28" "d_31" "d_33" "d_36" "d_39" "d_4" 
## [11] "d_7"  "d_11" "d_12" "d_19" "d_23" "d_26" "d_29" "d_32" "d_35" "d_38"
## [21] "d_40" "d_8"  "d_42" "d_44" "d_46" "d_50" "d_52" "d_54" "d_56" "d_58"
## [31] "d_61" "d_63" "d_65" "d_68" "d_15" "d_41" "d_43" "d_45" "d_47" "d_51"
## [41] "d_53" "d_67" "d_49" "d_14" "d_48" "d_37" "d_55" "d_57" "d_64" "d_66"
## [51] "d_16" "d_2"  "d_5"  "d_60" "d_62" "d_34" "d_10" "d_13" "d_21" "d_3" 
## [61] "d_6"  "d_17" "d_20"

Ind.cv.B <- createCV(denver.model.B, groups = 10, #min.dist = .1,
                     subset = site_idsB)

ID.cv.B <- sapply(split(denver.model.B$obs$ID, Ind.cv.B), unique)
print(sapply(ID.cv.B, length))

##  0  1  2  3  4  5  6  7  8  9 10 
##  1  7  7  7  6  6  6  6  6  6  6

table(Ind.cv.B)

## Ind.cv.B
##   0   1   2   3   4   5   6   7   8   9  10 
## 231  66  76  82  67  57  85  72  78  54  75

I.col.B <- apply(sapply(ID.cv.B, function(x) denver.model.B$locations$ID%in% x), 1,
                 function(x) if(sum(x)==1) which(x) else 0)
names(I.col.B) <- denver.model.B$locations$ID
print(I.col.B)

## central     d_1    d_10    d_11    d_12    d_13    d_14    d_15    d_16    d_17 
##       1       6       2       3       4       4       4       9       4      10 
##    d_18    d_19     d_2    d_20    d_21    d_22    d_23    d_25    d_26    d_28 
##       8       5      11       3       2      11       4       9      11       9 
##    d_29     d_3    d_31    d_32    d_33    d_34    d_35    d_36    d_37    d_38 
##       2       4       6       8       5       5      10       8       3       8 
##    d_39     d_4    d_40    d_41    d_42    d_43    d_44    d_45    d_46    d_47 
##       6       7       7       3       7       8      11       3       9       5 
##    d_48    d_49     d_5    d_50    d_51    d_52    d_53    d_54    d_55    d_56 
##       9       4       6       9      10       3       7       5       7       6 
##    d_57    d_58     d_6    d_60    d_61    d_62    d_63    d_64    d_65    d_66 
##       2       2      11       5       2      11      10       3       8       6 
##    d_67    d_68     d_7     d_8    d_24    d_27    d_30    d_59     d_9 
##       7      10       2      10       0       0       0       0       0

par(mfrow=c(1,1))
plot(denver.model.B$locations$long,
     denver.model.B$locations$lat,
     pch=23+floor(I.col.B/max(I.col.B)+.5), bg=I.col.B,
     xlab="Longitude", ylab="Latitude")

map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)

## [[1]]
## NULL

ID starting conditions using previous model without CV:

x.init.B.cv <- coef(est.denver.model.B, pars="cov")[,c("par","init")]
x.init.B.cv

Run the model with cross validation

est.denver.B.cv <- estimateCV(denver.model.B, x.init.B.cv, Ind.cv.B)

## 
## ***************************

## Estimation of CV-set 1/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1325.1  |proj g|=       14.284
## At iterate    10  f =      -1325.5  |proj g|=      0.089734
## At iterate    20  f =      -1325.6  |proj g|=      0.071136
## 
## iterations 25
## function evaluations 41
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0380497
## final function value -1325.55
## 
## F = -1325.55
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1325.550743 
## converged

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1020.4  |proj g|=           20
## At iterate    10  f =      -1323.1  |proj g|=        3.4163
## At iterate    20  f =      -1325.5  |proj g|=       0.22009
## At iterate    30  f =      -1325.5  |proj g|=       0.21144
## At iterate    40  f =      -1325.6  |proj g|=      0.028442
## 
## iterations 41
## function evaluations 50
## segments explored during Cauchy searches 46
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0159431
## final function value -1325.55
## 
## F = -1325.55
## final  value -1325.550741 
## converged

## 

## 
## ***************************

## Estimation of CV-set 2/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1337.4  |proj g|=       10.182
## At iterate    10  f =        -1338  |proj g|=       0.20161
## At iterate    20  f =        -1338  |proj g|=       0.50961
## At iterate    30  f =      -1338.1  |proj g|=       0.10069
## At iterate    40  f =      -1338.1  |proj g|=       0.17518
## At iterate    50  f =      -1338.1  |proj g|=       0.27559
## At iterate    60  f =      -1338.1  |proj g|=       0.06003
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## 
## iterations 69
## function evaluations 111
## segments explored during Cauchy searches 72
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00858394
## final function value -1338.13
## 
## F = -1338.13
## final  value -1338.131602 
## converged

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1032.7  |proj g|=           20
## At iterate    10  f =      -1335.3  |proj g|=        4.0886
## At iterate    20  f =      -1338.1  |proj g|=       0.59767
## At iterate    30  f =      -1338.1  |proj g|=        1.0666
## At iterate    40  f =      -1338.1  |proj g|=      0.045435
## At iterate    50  f =      -1338.1  |proj g|=      0.078825
## 
## iterations 50
## function evaluations 66
## segments explored during Cauchy searches 55
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0788252
## final function value -1338.13
## 
## F = -1338.13
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1338.131578 
## converged

## 

## 
## ***************************

## Estimation of CV-set 3/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1322.9  |proj g|=       9.3055
## At iterate    10  f =      -1323.3  |proj g|=       0.44233
## At iterate    20  f =      -1323.4  |proj g|=      0.077529
## 
## iterations 27
## function evaluations 46
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0242257
## final function value -1323.41
## 
## F = -1323.41
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1323.412118 
## converged

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1025.7  |proj g|=           20
## At iterate    10  f =      -1319.3  |proj g|=        20.067
## At iterate    20  f =      -1322.9  |proj g|=       0.17057
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## 
## iterations 29
## function evaluations 74
## segments explored during Cauchy searches 35
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0414085
## final function value -1323.41
## 
## F = -1323.41
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1323.412120 
## converged

## 

## 
## ***************************

## Estimation of CV-set 4/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1321.3  |proj g|=       16.868
## At iterate    10  f =      -1322.2  |proj g|=       0.14931
## At iterate    20  f =      -1322.2  |proj g|=       0.48659
## At iterate    30  f =      -1322.2  |proj g|=      0.019047
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1322.239329 
## stopped after 32 iterations

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1019.9  |proj g|=           20
## At iterate    10  f =      -1319.5  |proj g|=        3.5891
## At iterate    20  f =      -1322.2  |proj g|=      0.033125
## 
## iterations 22
## function evaluations 29
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0418117
## final function value -1322.23
## 
## F = -1322.23
## final  value -1322.233610 
## converged

## 

## 
## ***************************

## Estimation of CV-set 5/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1346.7  |proj g|=        10.31
## At iterate    10  f =      -1347.3  |proj g|=       0.24337
## At iterate    20  f =      -1347.4  |proj g|=       0.78481
## At iterate    30  f =      -1347.5  |proj g|=       0.48791
## At iterate    40  f =      -1347.5  |proj g|=      0.076001
## At iterate    50  f =      -1347.5  |proj g|=      0.053906
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## 
## iterations 53
## function evaluations 93
## segments explored during Cauchy searches 54
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0640349
## final function value -1347.52
## 
## F = -1347.52
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1347.517035 
## converged

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1043.4  |proj g|=           20
## At iterate    10  f =      -1345.7  |proj g|=         3.439
## At iterate    20  f =      -1347.5  |proj g|=        0.2081
## At iterate    30  f =      -1347.5  |proj g|=        1.2576
## At iterate    40  f =      -1347.5  |proj g|=       0.13001
## 
## iterations 46
## function evaluations 63
## segments explored during Cauchy searches 51
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00663204
## final function value -1347.52
## 
## F = -1347.52
## final  value -1347.517059 
## converged

## 

## 
## ***************************

## Estimation of CV-set 6/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1320.9  |proj g|=       10.182
## At iterate    10  f =      -1321.2  |proj g|=       0.27646
## At iterate    20  f =      -1321.2  |proj g|=      0.066922
## At iterate    30  f =      -1321.2  |proj g|=      0.035871
## 
## iterations 30
## function evaluations 36
## segments explored during Cauchy searches 31
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0358709
## final function value -1321.18
## 
## F = -1321.18
## final  value -1321.183539 
## converged

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1008.5  |proj g|=           20
## At iterate    10  f =      -1319.2  |proj g|=        3.6121
## At iterate    20  f =      -1321.2  |proj g|=      0.097296
## 
## iterations 24
## function evaluations 30
## segments explored during Cauchy searches 29
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0223154
## final function value -1321.17
## 
## F = -1321.17
## final  value -1321.171547 
## converged

## 

## 
## ***************************

## Estimation of CV-set 7/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1330.9  |proj g|=       2.1047
## At iterate    10  f =        -1332  |proj g|=        1.8288
## At iterate    20  f =      -1332.1  |proj g|=       0.66345
## At iterate    30  f =      -1332.3  |proj g|=       0.14817
## 
## iterations 38
## function evaluations 43
## segments explored during Cauchy searches 39
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 2
## norm of the final projected gradient 0.000450916
## final function value -1332.32
## 
## F = -1332.32
## final  value -1332.320536 
## converged

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1017.7  |proj g|=           20
## At iterate    10  f =      -1328.6  |proj g|=        2.9518
## At iterate    20  f =      -1332.2  |proj g|=       0.43645
## At iterate    30  f =      -1332.3  |proj g|=       0.45755
## At iterate    40  f =      -1332.3  |proj g|=      0.021828
## 
## iterations 43
## function evaluations 58
## segments explored during Cauchy searches 48
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0851882
## final function value -1332.32
## 
## F = -1332.32
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1332.320410 
## converged

## 

## 
## ***************************

## Estimation of CV-set 8/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1308.2  |proj g|=        12.34
## At iterate    10  f =      -1309.5  |proj g|=       0.24255
## At iterate    20  f =      -1309.6  |proj g|=      0.047146
## At iterate    30  f =      -1309.6  |proj g|=      0.080695
## At iterate    40  f =      -1309.6  |proj g|=      0.011327
## 
## iterations 40
## function evaluations 58
## segments explored during Cauchy searches 40
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0113275
## final function value -1309.59
## 
## F = -1309.59
## final  value -1309.594791 
## converged

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1003.1  |proj g|=           20
## At iterate    10  f =      -1307.2  |proj g|=        2.8981
## At iterate    20  f =      -1309.3  |proj g|=        1.2754
## At iterate    30  f =      -1309.4  |proj g|=        1.6429
## At iterate    40  f =      -1309.5  |proj g|=       0.22026
## At iterate    50  f =      -1309.5  |proj g|=        1.8767
## At iterate    60  f =      -1309.6  |proj g|=      0.094154
## At iterate    70  f =      -1309.6  |proj g|=      0.012398
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## 
## iterations 72
## function evaluations 122
## segments explored during Cauchy searches 78
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0509512
## final function value -1309.59
## 
## F = -1309.59
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1309.594828 
## converged

## 

## 
## ***************************

## Estimation of CV-set 9/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=        -1353  |proj g|=         8.74
## At iterate    10  f =      -1353.2  |proj g|=       0.21479
## 
## iterations 13
## function evaluations 16
## segments explored during Cauchy searches 13
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0334872
## final function value -1353.24
## 
## F = -1353.24
## final  value -1353.238180 
## converged

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1046.2  |proj g|=           20
## At iterate    10  f =      -1350.5  |proj g|=        3.0881
## At iterate    20  f =      -1353.2  |proj g|=       0.13529
## 
## iterations 23
## function evaluations 30
## segments explored during Cauchy searches 28
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0461884
## final function value -1353.2
## 
## F = -1353.2
## final  value -1353.203269 
## converged

## 

## 
## ***************************

## Estimation of CV-set 10/10

## Optimisation using starting value 1/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1334.7  |proj g|=       9.2237
## At iterate    10  f =      -1335.3  |proj g|=       0.31064
## At iterate    20  f =      -1335.3  |proj g|=        0.1466
## 
## iterations 25
## function evaluations 38
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0154309
## final function value -1335.29
## 
## F = -1335.29
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1335.293048 
## converged

## 

## Optimisation using starting value 2/2

## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1021.3  |proj g|=           20
## At iterate    10  f =      -1333.1  |proj g|=        2.4865
## At iterate    20  f =      -1335.1  |proj g|=       0.69007
## At iterate    30  f =      -1335.2  |proj g|=       0.98087
## At iterate    40  f =      -1335.3  |proj g|=      0.024394
## 
## iterations 40
## function evaluations 49
## segments explored during Cauchy searches 45
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0243941
## final function value -1335.29
## 
## F = -1335.29
## final  value -1335.292899 
## converged

## 

print(est.denver.B.cv)

## Cross-validation parameter estimation for STmodel
##   with 10 CV-groups and 2 starting points.
##   Results: 10 converged, 0 not converged.
## 
## No fixed parameters.
## 
## Estimated function values and convergence info:
##       value convergence conv    eigen.min eigen.all.min
## 1  1325.551        TRUE TRUE 0.3083541248            NA
## 2  1338.132        TRUE TRUE 0.0013147477            NA
## 3  1323.412        TRUE TRUE 1.0261104123            NA
## 4  1322.234        TRUE TRUE 0.0019628702            NA
## 5  1347.517        TRUE TRUE 0.0014229949            NA
## 6  1321.184        TRUE TRUE 0.0903619032            NA
## 7  1332.321        TRUE TRUE 0.0008020935            NA
## 8  1309.595        TRUE TRUE 0.1956246001            NA
## 9  1353.238        TRUE TRUE 0.2799606817            NA
## 10 1335.293        TRUE TRUE 0.5380999282            NA

par(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.B, pars="all"),
     plotCI((1:length(par))+.3, par, uiw=1.96*sd,
            col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.B.cv, "all", boxwex=.4, col="grey", add=TRUE)

#' Save the results as an .rdata object
save(denver.data.B, denver.model.B, est.denver.model.B, est.denver.B.cv,
     file = here::here("Results", "Denver_ST_Model_B.rdata"))

3.2.4 Prediction using the CV model

Making predictions using the CV model. Printing out the CV summary statistics as well

pred.B.cv <- predictCV(denver.model.B, est.denver.B.cv, LTA = T)
pred.B.cv.log <- predictCV(denver.model.B, est.denver.B.cv, LTA = T,
                           transform="unbiased")


head(pred.B.cv$pred.all$EX)

##            central       d_1      d_10      d_11      d_12      d_13      d_14
## 2008-12-29      NA 0.2256719 0.2715704 0.4040500 0.3325563 0.3381283 0.3552283
## 2009-01-05      NA 0.2454387 0.2867163 0.4125644 0.3424156 0.3470399 0.3569286
## 2009-01-12      NA 0.3202412 0.3522732 0.4671434 0.3973428 0.4010627 0.4036430
## 2009-01-19      NA 0.2611628 0.2956337 0.4101100 0.3424710 0.3453806 0.3404619
## 2009-01-26      NA 0.2625262 0.2945710 0.4039017 0.3370371 0.3392808 0.3265794
## 2009-02-02      NA 0.2774678 0.3063329 0.4094768 0.3426068 0.3443714 0.3235264
##                  d_15      d_16      d_17      d_18      d_19       d_2
## 2008-12-29 0.08474891 0.4246954 0.2028410 0.3308776 0.2007175 0.3439717
## 2009-01-05 0.12050759 0.4430645 0.2272220 0.3412991 0.2302263 0.3775853
## 2009-01-12 0.21196409 0.5071419 0.3034162 0.4055201 0.3213495 0.4671328
## 2009-01-19 0.16712856 0.4627590 0.2527525 0.3395599 0.2637587 0.4242432
## 2009-01-26 0.18303365 0.4698510 0.2585887 0.3330385 0.2743600 0.4444420
## 2009-02-02 0.21192172 0.4905654 0.2764793 0.3404257 0.3011394 0.4815097
##                 d_20      d_21      d_22      d_23      d_25      d_26
## 2008-12-29 0.3708507 0.2226343 0.2735033 0.3692873 0.2239682 0.2836085
## 2009-01-05 0.3803355 0.2426663 0.2926294 0.3848083 0.2417591 0.3000840
## 2009-01-12 0.4359096 0.3129022 0.3667768 0.4456929 0.3156517 0.3718684
## 2009-01-19 0.3799220 0.2604923 0.3065849 0.3973926 0.2541472 0.3099290
## 2009-01-26 0.3748357 0.2629673 0.3065892 0.3994628 0.2547692 0.3091289
## 2009-02-02 0.3816306 0.2773724 0.3197456 0.4137345 0.2701712 0.3226977
##                 d_28      d_29       d_3      d_31      d_32      d_33
## 2008-12-29 0.2037387 0.2646549 0.3040163 0.3402394 0.2963519 0.1975749
## 2009-01-05 0.2205081 0.2741179 0.3174851 0.3546011 0.3047909 0.2182837
## 2009-01-12 0.2933147 0.3338391 0.3759481 0.4245347 0.3669623 0.3005877
## 2009-01-19 0.2305899 0.2710521 0.3244478 0.3617353 0.2988150 0.2341560
## 2009-01-26 0.2297874 0.2633712 0.3221305 0.3611373 0.2898981 0.2358924
## 2009-02-02 0.2435023 0.2679122 0.3304864 0.3763881 0.2946226 0.2537815
##                 d_34      d_35      d_36      d_37      d_38      d_39
## 2008-12-29 0.3537773 0.3138491 0.3399919 0.3601579 0.2811222 0.3474737
## 2009-01-05 0.3678469 0.3254832 0.3488600 0.3679288 0.2988888 0.3520544
## 2009-01-12 0.4438244 0.3891428 0.4115710 0.4215859 0.3705696 0.4120839
## 2009-01-19 0.3717432 0.3264221 0.3441967 0.3632544 0.3122969 0.3391408
## 2009-01-26 0.3688713 0.3209440 0.3364089 0.3551737 0.3138032 0.3280427
## 2009-02-02 0.3834977 0.3284855 0.3427212 0.3581371 0.3296518 0.3323407
##                  d_4      d_40      d_41      d_42      d_43      d_44
## 2008-12-29 0.2382955 0.2492708 0.3870883 0.2337331 0.3137349 0.3119008
## 2009-01-05 0.2686225 0.2713930 0.3883383 0.2593849 0.3231533 0.3263064
## 2009-01-12 0.3808060 0.3754023 0.4356365 0.3668958 0.3865030 0.3960576
## 2009-01-19 0.2981158 0.2846203 0.3713036 0.2795488 0.3199534 0.3321672
## 2009-01-26 0.3073177 0.2858647 0.3577746 0.2841226 0.3132742 0.3295439
## 2009-02-02 0.3403762 0.3111443 0.3560063 0.3125931 0.3210603 0.3414555
##                 d_45      d_46      d_47      d_48      d_49       d_5
## 2008-12-29 0.3965320 0.2567529 0.1750927 0.3120921 0.3362977 0.2812456
## 2009-01-05 0.3963727 0.2742614 0.2026760 0.3250762 0.3459776 0.2983658
## 2009-01-12 0.4424060 0.3478701 0.2918032 0.3941431 0.4008662 0.3706735
## 2009-01-19 0.3771173 0.2860787 0.2320709 0.3277827 0.3462545 0.3094279
## 2009-01-26 0.3631061 0.2864096 0.2403103 0.3235076 0.3415375 0.3091272
## 2009-02-02 0.3614669 0.3015152 0.2644451 0.3339626 0.3484125 0.3230546
##                 d_50      d_51      d_52      d_53      d_54      d_55
## 2008-12-29 0.2118696 0.2859845 0.3952992 0.1619867 0.2032670 0.2547967
## 2009-01-05 0.2300339 0.2994959 0.4021042 0.1961947 0.2270731 0.2729522
## 2009-01-12 0.3043488 0.3650713 0.4549500 0.3121661 0.3125292 0.3728944
## 2009-01-19 0.2433689 0.3043440 0.3961372 0.2330581 0.2493588 0.2778413
## 2009-01-26 0.2446715 0.3009760 0.3880805 0.2455238 0.2545202 0.2745069
## 2009-02-02 0.2609549 0.3107765 0.3917193 0.2814335 0.2760431 0.2948136
##                 d_56      d_57      d_58       d_6      d_60      d_61
## 2008-12-29 0.2854202 0.4341711 0.2565010 0.2954464 0.2226554 0.2857551
## 2009-01-05 0.3015254 0.4326568 0.2746232 0.3090520 0.2443350 0.3015237
## 2009-01-12 0.3728553 0.4816223 0.3430821 0.3779759 0.3276838 0.3675802
## 2009-01-19 0.3107126 0.4085735 0.2891807 0.3132019 0.2624518 0.3111780
## 2009-01-26 0.3096394 0.3913963 0.2906024 0.3096097 0.2656235 0.3099508
## 2009-02-02 0.3229555 0.3874346 0.3045200 0.3204431 0.2852508 0.3210304
##                 d_62      d_63      d_64      d_65      d_66      d_67
## 2008-12-29 0.3313692 0.2867062 0.3873173 0.2832015 0.2454683 0.2485100
## 2009-01-05 0.3458785 0.3005263 0.3861333 0.3000520 0.2633762 0.2763734
## 2009-01-12 0.4156277 0.3663932 0.4308383 0.3707025 0.3364746 0.3862913
## 2009-01-19 0.3515107 0.3059206 0.3635800 0.3111593 0.2760240 0.3017601
## 2009-01-26 0.3483175 0.3027505 0.3466204 0.3110289 0.2765243 0.3097724
## 2009-02-02 0.3592168 0.3126757 0.3407722 0.3247696 0.2912599 0.3424822
##                 d_68       d_7       d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.2816600 0.2675658 0.2719639   NA   NA   NA   NA  NA
## 2009-01-05 0.2998819 0.2833328 0.2883306   NA   NA   NA   NA  NA
## 2009-01-12 0.3700401 0.3493903 0.3567725   NA   NA   NA   NA  NA
## 2009-01-19 0.3136151 0.2929949 0.2989293   NA   NA   NA   NA  NA
## 2009-01-26 0.3141158 0.2917834 0.2984695   NA   NA   NA   NA  NA
## 2009-02-02 0.3272234 0.3028901 0.3112074   NA   NA   NA   NA  NA

tail(pred.B.cv$pred.all$EX)

##            central       d_1      d_10     d_11      d_12      d_13      d_14
## 2020-11-02      NA 1.0592974 1.0228608 1.061747 1.0395015 1.0352110 0.9430054
## 2020-11-09      NA 0.6743546 0.6436968 0.670476 0.6555533 0.6414596 0.5832082
## 2020-11-16      NA 1.0748031 1.0491068 1.086222 1.0686527 1.0601094 1.0035757
## 2020-11-23      NA 1.0128161 0.9815864 1.026519 1.0057695 0.9995560 0.9514267
## 2020-11-30      NA 1.2072519 1.1746121 1.236153 1.2102909 1.2097030 1.1687079
## 2020-12-07      NA        NA        NA       NA        NA        NA        NA
##                 d_15      d_16      d_17      d_18      d_19       d_2
## 2020-11-02 1.1419553 1.1275054 1.1310229 1.0238101 1.1039037 1.2837506
## 2020-11-09 0.7476813 0.7102194 0.7356171 0.6491306 0.7212576 0.8284186
## 2020-11-16 1.1288283 1.1099571 1.1315912 1.0667230 1.1094668 1.2016631
## 2020-11-23 1.0364911 1.0045859 1.0578360 1.0122618 1.0301734 1.0992655
## 2020-11-30 1.1975460 1.1997819 1.2452609 1.2195107 1.1992809 1.2551670
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_20      d_21      d_22      d_23      d_25     d_26      d_28
## 2020-11-02 1.0191267 1.0805214 1.1344515 1.0770611 0.9911015 1.063403 0.9644699
## 2020-11-09 0.6419187 0.6906482 0.7454757 0.6860287 0.6114598 0.662137 0.5896073
## 2020-11-16 1.0581193 1.0906546 1.1559272 1.0892098 1.0168345 1.078312 0.9941662
## 2020-11-23 0.9940127 1.0211557 1.0949061 1.0059062 0.9483883 1.018356 0.9300305
## 2020-11-30 1.2055783 1.2123088 1.2955459 1.1892691 1.1432666 1.222471 1.1271004
## 2020-12-07        NA        NA        NA        NA        NA       NA        NA
##                 d_29       d_3      d_31      d_32      d_33      d_34
## 2020-11-02 0.9756170 1.0777925 0.9988789 0.9793660 1.0054510 1.0214732
## 2020-11-09 0.6090801 0.6853581 0.6207785 0.6114122 0.6246790 0.6480789
## 2020-11-16 1.0240521 1.0964536 1.0265601 1.0334813 1.0267347 1.0591085
## 2020-11-23 0.9700332 1.0282825 0.9582950 0.9824494 0.9663557 1.0026145
## 2020-11-30 1.1776068 1.2287002 1.1473745 1.1940734 1.1596187 1.2034563
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_35      d_36      d_37      d_38      d_39       d_4
## 2020-11-02 1.0415438 1.0048275 1.0566150 1.0367478 0.9030595 1.0685735
## 2020-11-09 0.6667445 0.6333597 0.6803382 0.6496136 0.5461058 0.6657496
## 2020-11-16 1.0822076 1.0529320 1.1005811 1.0523396 0.9632485 1.0695919
## 2020-11-23 1.0240669 0.9997173 1.0412693 0.9802670 0.9176003 1.0111316
## 2020-11-30 1.2313727 1.2084478 1.2506371 1.1743400 1.1272262 1.2000965
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_40      d_41      d_42      d_43      d_44      d_45
## 2020-11-02 0.8413985 0.9707494 1.0133947 0.9667202 1.0469778 0.9281184
## 2020-11-09 0.4871106 0.6037283 0.6257857 0.5950959 0.6599757 0.5622223
## 2020-11-16 0.8729626 1.0334131 1.0379072 1.0126210 1.0738734 0.9901673
## 2020-11-23 0.8106980 0.9794114 0.9843721 0.9539198 1.0150174 0.9377863
## 2020-11-30 0.9991102 1.1997688 1.1778619 1.1613094 1.2177293 1.1555613
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_46      d_47      d_48      d_49       d_5      d_50
## 2020-11-02 1.0114037 1.0703124 1.0200884 1.0144921 1.0151890 0.9595869
## 2020-11-09 0.6285872 0.6935142 0.6472952 0.6318189 0.6376615 0.5799605
## 2020-11-16 1.0296364 1.0842739 1.0581806 1.0464654 1.0321577 0.9765684
## 2020-11-23 0.9687930 1.0127037 1.0034845 0.9787747 0.9695431 0.9127071
## 2020-11-30 1.1631952 1.1797399 1.2050546 1.1762742 1.1639134 1.1030276
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_51      d_52      d_53      d_54      d_55      d_56
## 2020-11-02 1.0180081 1.0472624 1.1533852 1.0267021 0.9534371 1.0170177
## 2020-11-09 0.6499887 0.6646837 0.7306014 0.6454562 0.5876489 0.6347523
## 2020-11-16 1.0606907 1.0842299 1.1325974 1.0453139 1.0116455 1.0350990
## 2020-11-23 0.9989803 1.0264275 1.0687857 0.9792002 0.9695331 0.9763999
## 2020-11-30 1.1899641 1.2363109 1.2482694 1.1734064 1.1756000 1.1742615
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_57      d_58       d_6      d_60     d_61      d_62      d_63
## 2020-11-02 0.9611157 1.0753301 1.0284772 1.0142264 1.053816 1.0684598 1.0344341
## 2020-11-09 0.6116437 0.6854853 0.6472296 0.6319410 0.667495 0.6844248 0.6548452
## 2020-11-16 1.0404840 1.0871688 1.0638627 1.0287715 1.069091 1.0941749 1.0645456
## 2020-11-23 1.0018217 1.0186686 1.0041998 0.9629806 1.007879 1.0384102 1.0046810
## 2020-11-30 1.2233153 1.2082711 1.2071578 1.1537165 1.207289 1.2432936 1.2054039
## 2020-12-07        NA        NA        NA        NA       NA        NA        NA
##                 d_64      d_65      d_66      d_67      d_68       d_7
## 2020-11-02 0.9675467 1.0558200 1.0256708 1.0136834 1.1363678 1.0344538
## 2020-11-09 0.6091437 0.6490888 0.6429818 0.6305675 0.7437944 0.6484312
## 2020-11-16 1.0417461 1.0617888 1.0468392 1.0376018 1.1498478 1.0503399
## 2020-11-23 0.9985844 1.0027639 0.9801774 0.9728781 1.0854197 0.9874180
## 2020-11-30 1.2229546 1.2085839 1.1764166 1.1494601 1.2794957 1.1866676
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                  d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.0464963   NA   NA   NA   NA  NA
## 2020-11-09 0.6604396   NA   NA   NA   NA  NA
## 2020-11-16 1.0660222   NA   NA   NA   NA  NA
## 2020-11-23 1.0015201   NA   NA   NA   NA  NA
## 2020-11-30 1.1976984   NA   NA   NA   NA  NA
## 2020-12-07        NA   NA   NA   NA   NA  NA

head(pred.B.cv.log$pred.all$EX)

##            central      d_1     d_10     d_11     d_12     d_13     d_14
## 2008-12-29      NA 1.285715 1.346294 1.537529 1.432267 1.440270 1.465110
## 2009-01-05      NA 1.311077 1.366507 1.550255 1.446051 1.452754 1.467191
## 2009-01-12      NA 1.412620 1.458780 1.636827 1.527324 1.533016 1.536977
## 2009-01-19      NA 1.331348 1.378201 1.545771 1.445474 1.449686 1.442573
## 2009-01-26      NA 1.332971 1.376532 1.535948 1.437400 1.440628 1.422446
## 2009-02-02      NA 1.352884 1.392658 1.544336 1.445244 1.447797 1.417930
##                d_15     d_16     d_17     d_18     d_19      d_2     d_20
## 2008-12-29 1.116118 1.570506 1.256995 1.428666 1.253541 1.448685 1.487322
## 2009-01-05 1.156607 1.599171 1.287688 1.443225 1.290773 1.497763 1.501089
## 2009-01-12 1.267234 1.704576 1.389325 1.538570 1.413621 1.637645 1.586493
## 2009-01-19 1.211556 1.630237 1.320438 1.440057 1.334269 1.568547 1.499804
## 2009-01-26 1.230881 1.641564 1.327962 1.430453 1.348290 1.600270 1.491947
## 2009-02-02 1.266877 1.675710 1.351773 1.440872 1.384725 1.660480 1.501925
##                d_21     d_22     d_23     d_25     d_26     d_28     d_29
## 2008-12-29 1.281998 1.350112 1.485854 1.282840 1.363825 1.257149 1.337016
## 2009-01-05 1.307619 1.375774 1.508671 1.305704 1.386068 1.278249 1.349400
## 2009-01-12 1.402463 1.481275 1.602984 1.405684 1.488836 1.374634 1.432135
## 2009-01-19 1.330610 1.394438 1.527082 1.321707 1.399109 1.290935 1.344736
## 2009-01-26 1.333709 1.394198 1.529990 1.322424 1.397744 1.289796 1.334248
## 2009-02-02 1.352904 1.412474 1.551785 1.342863 1.416650 1.307524 1.340166
##                 d_3     d_31     d_32     d_33     d_34     d_35     d_36
## 2008-12-29 1.391968 1.441786 1.380182 1.249607 1.460869 1.404571 1.441747
## 2009-01-05 1.410446 1.462301 1.391486 1.275450 1.481214 1.420644 1.454179
## 2009-01-12 1.494994 1.567904 1.480375 1.384574 1.597802 1.513681 1.547908
## 2009-01-19 1.419655 1.472210 1.382561 1.295350 1.486416 1.421387 1.446750
## 2009-01-26 1.416132 1.471117 1.370055 1.297409 1.481935 1.413403 1.435282
## 2009-02-02 1.427833 1.493555 1.376364 1.320676 1.503596 1.423934 1.444183
##                d_37     d_38     d_39      d_4     d_40     d_41     d_42
## 2008-12-29 1.471503 1.359322 1.452254 1.302283 1.316655 1.511670 1.296355
## 2009-01-05 1.482580 1.383297 1.458582 1.342038 1.345762 1.513150 1.329698
## 2009-01-12 1.563930 1.485725 1.548503 1.501024 1.492935 1.586060 1.480289
## 2009-01-19 1.475013 1.401327 1.439319 1.381632 1.363111 1.486934 1.356216
## 2009-01-26 1.462899 1.403201 1.423227 1.394185 1.364594 1.466708 1.362219
## 2009-02-02 1.467051 1.425431 1.429195 1.440870 1.399360 1.463928 1.401389
##                d_43     d_44     d_45     d_46     d_47     d_48     d_49
## 2008-12-29 1.404384 1.402961 1.526013 1.325594 1.221827 1.401019 1.437636
## 2009-01-05 1.417273 1.422895 1.525356 1.348839 1.255698 1.419152 1.451211
## 2009-01-12 1.509587 1.525289 1.596833 1.451711 1.372464 1.520464 1.532715
## 2009-01-19 1.412098 1.430571 1.495604 1.364592 1.292652 1.422704 1.450953
## 2009-01-26 1.402459 1.426572 1.474549 1.364934 1.303154 1.416522 1.443883
## 2009-02-02 1.413237 1.443474 1.471944 1.385621 1.334834 1.431318 1.453660
##                 d_5     d_50     d_51     d_52     d_53     d_54     d_55
## 2008-12-29 1.359190 1.267413 1.365974 1.524133 1.206605 1.256740 1.323951
## 2009-01-05 1.382338 1.290483 1.384200 1.534124 1.248274 1.286710 1.347862
## 2009-01-12 1.485689 1.389885 1.477680 1.616990 1.401450 1.401207 1.489195
## 2009-01-19 1.397182 1.307538 1.390349 1.524322 1.294608 1.315193 1.353902
## 2009-01-26 1.396559 1.309137 1.385460 1.511839 1.310641 1.321804 1.349183
## 2009-02-02 1.415985 1.330544 1.398939 1.517154 1.358396 1.350406 1.376693
##                d_56     d_57     d_58      d_6     d_60     d_61     d_62
## 2008-12-29 1.364876 1.584005 1.326158 1.380065 1.281344 1.365527 1.430542
## 2009-01-05 1.386712 1.581223 1.350082 1.398554 1.309114 1.386892 1.451018
## 2009-01-12 1.488934 1.660220 1.445434 1.497957 1.422603 1.481282 1.555433
## 2009-01-19 1.398978 1.542985 1.369336 1.403696 1.332526 1.399791 1.458513
## 2009-01-26 1.397275 1.516481 1.371080 1.398416 1.336562 1.397867 1.453607
## 2009-02-02 1.415845 1.510311 1.390135 1.413460 1.362897 1.413278 1.469341
##                d_63     d_64     d_65     d_66     d_67     d_68      d_7
## 2008-12-29 1.366960 1.512016 1.362151 1.311422 1.315654 1.360079 1.340914
## 2009-01-05 1.385627 1.509817 1.384907 1.334807 1.352481 1.384735 1.361892
## 2009-01-12 1.479634 1.578468 1.485922 1.435739 1.509280 1.485040 1.454581
## 2009-01-19 1.392543 1.475494 1.399734 1.351281 1.386676 1.403299 1.374569
## 2009-01-26 1.387921 1.450439 1.399313 1.351762 1.397612 1.403785 1.372700
## 2009-02-02 1.401599 1.441796 1.418489 1.371673 1.443908 1.422138 1.387871
##                 d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 1.346955   NA   NA   NA   NA  NA
## 2009-01-05 1.368831   NA   NA   NA   NA  NA
## 2009-01-12 1.465467   NA   NA   NA   NA  NA
## 2009-01-19 1.382841   NA   NA   NA   NA  NA
## 2009-01-26 1.381992   NA   NA   NA   NA  NA
## 2009-02-02 1.399542   NA   NA   NA   NA  NA

tail(pred.B.cv.log$pred.all$EX)

##            central      d_1     d_10     d_11     d_12     d_13     d_14
## 2020-11-02      NA 2.899503 2.794215 2.904608 2.842982 2.829757 2.581567
## 2020-11-09      NA 1.972851 1.912181 1.963776 1.936143 1.908337 1.801091
## 2020-11-16      NA 2.944390 2.867957 2.975923 2.926122 2.900149 2.741870
## 2020-11-23      NA 2.767585 2.680783 2.803544 2.747680 2.729643 2.602450
## 2020-11-30      NA 3.362062 3.251912 3.457825 3.371360 3.368125 3.234170
## 2020-12-07      NA       NA       NA       NA       NA       NA       NA
##                d_15     d_16     d_17     d_18     d_19      d_2     d_20
## 2020-11-02 3.145770 3.101308 3.114198 2.796847 3.029272 3.627302 2.782459
## 2020-11-09 2.120700 2.042818 2.096767 1.922487 2.065863 2.300051 1.907834
## 2020-11-16 3.104689 3.046361 3.115183 2.918621 3.045665 3.340231 2.892463
## 2020-11-23 2.831069 2.741593 2.893688 2.763885 2.813628 3.014971 2.712944
## 2020-11-30 3.326239 3.332671 3.490455 3.400611 3.332469 3.523719 3.352551
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_21     d_22     d_23     d_25     d_26     d_28     d_29
## 2020-11-02 2.960609 3.126655 2.949135 2.706198 2.910953 2.635964 2.666323
## 2020-11-09 2.004469 2.118615 1.994257 1.851257 1.948362 1.811850 1.847847
## 2020-11-16 2.990173 3.193375 2.984204 2.776691 2.953610 2.715368 2.798094
## 2020-11-23 2.789498 3.004181 2.745579 2.593196 2.781582 2.546880 2.651032
## 2020-11-30 3.377457 3.671754 3.298255 3.151575 3.411527 3.102075 3.262947
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_3     d_31     d_32     d_33     d_34     d_35     d_36
## 2020-11-02 2.952940 2.727177 2.675206 2.746483 2.789847 2.846940 2.744197
## 2020-11-09 1.994034 1.868337 1.851283 1.876524 1.920267 1.956733 1.892365
## 2020-11-16 3.007579 2.803318 2.823133 2.805100 2.896363 2.964346 2.878585
## 2020-11-23 2.809275 2.618494 2.682643 2.640871 2.737406 2.796909 2.729371
## 2020-11-30 3.432821 3.163955 3.315126 3.204329 3.346730 3.441461 3.363126
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_37     d_38     d_39      d_4     d_40     d_41     d_42
## 2020-11-02 2.889273 2.832845 2.480460 2.926076 2.330771 2.651326 2.768458
## 2020-11-09 1.982918 1.923129 1.735635 1.955587 1.635197 1.836534 1.878613
## 2020-11-16 3.018475 2.876513 2.633963 2.928462 2.405023 2.822165 2.836581
## 2020-11-23 2.844743 2.676456 2.516583 2.762262 2.259916 2.673893 2.688806
## 2020-11-30 3.507704 3.249937 3.103936 3.337194 2.728784 3.333471 3.263182
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_43     d_44     d_45     d_46     d_47     d_48     d_49
## 2020-11-02 2.641400 2.864568 2.541062 2.763561 2.929324 2.786918 2.770918
## 2020-11-09 1.821192 1.944860 1.762137 1.884506 2.009418 1.919577 1.889473
## 2020-11-16 2.764654 2.941595 2.703134 2.814360 2.970015 2.895074 2.860008
## 2020-11-23 2.607003 2.773316 2.565269 2.648435 2.765014 2.741189 2.672719
## 2020-11-30 3.208040 3.396618 3.189805 3.217173 3.268114 3.353800 3.256439
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_5     d_50     d_51     d_52     d_53     d_54     d_55
## 2020-11-02 2.774650 2.624441 2.780413 2.862679 3.182716 2.804348 2.607647
## 2020-11-09 1.901951 1.795357 1.924008 1.952325 2.085071 1.915152 1.808526
## 2020-11-16 2.821729 2.669341 2.900926 2.969835 3.116602 2.856557 2.763377
## 2020-11-23 2.650626 2.504396 2.727318 2.803130 2.924035 2.673937 2.649507
## 2020-11-30 3.219768 3.029810 3.301504 3.458176 3.499305 3.247513 3.256185
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_56     d_57     d_58      d_6     d_60     d_61     d_62
## 2020-11-02 2.779868 2.628354 2.944689 2.811351 2.769972 2.883928 2.928094
## 2020-11-09 1.896521 1.852884 1.993747 1.919744 1.889710 1.959500 1.993898
## 2020-11-16 2.830183 2.844904 2.979170 2.911560 2.810089 2.927741 3.003283
## 2020-11-23 2.668997 2.737092 2.782011 2.742786 2.631290 2.753983 2.840241
## 2020-11-30 3.253423 3.416094 3.363174 3.360054 3.184646 3.362102 3.486146
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_63     d_64     d_65     d_66     d_67     d_68      d_7
## 2020-11-02 2.827670 2.643929 2.887731 2.801725 2.767600 3.130032 2.827912
## 2020-11-09 1.934203 1.847262 1.922346 1.910622 1.886487 2.113405 1.922014
## 2020-11-16 2.913377 2.846945 2.904164 2.861252 2.834017 3.171711 2.872630
## 2020-11-23 2.744084 2.726770 2.737672 2.676897 2.656486 2.973804 2.697526
## 2020-11-30 3.354308 3.413058 3.363550 3.257765 3.169908 3.611032 3.292651
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.861798   NA   NA   NA   NA  NA
## 2020-11-09 1.944927   NA   NA   NA   NA  NA
## 2020-11-16 2.917491   NA   NA   NA   NA  NA
## 2020-11-23 2.735245   NA   NA   NA   NA  NA
## 2020-11-30 3.328343   NA   NA   NA   NA  NA
## 2020-12-07       NA   NA   NA   NA   NA  NA

summary(pred.B.cv)

## Cross-validation predictions for STmodel with 10 CV-groups.
##   Predictions for 712 observations.
##   Temporal averages for 63 locations.
## 
## RMSE:
##              EX.mu EX.mu.beta         EX
## obs     0.13672884 0.13672884 0.10865488
## average 0.08203204 0.08203204 0.06010949
## 
## R2:
##             EX.mu EX.mu.beta        EX
## obs     0.7318650  0.7318650 0.8306709
## average 0.7028679  0.7028679 0.8404603
## 
## Coverage of 95% prediction intervals:
##                EX
## obs     0.9550562
## average 0.9206349

summary(pred.B.cv.log)

## Cross-validation predictions for STmodel with 10 CV-groups.
##   Predictions for 712 observations.
##   Temporal averages for 63 locations.
## 
## RMSE:
##              EX.mu EX.mu.beta         EX    EX.pred
## obs     0.19696507 0.19696507 0.15055972 0.15053352
## average 0.09884076 0.09884076 0.06180058 0.06231628
## 
## R2:
##             EX.mu EX.mu.beta        EX   EX.pred
## obs     0.6975357  0.6975357 0.8232688 0.8233303
## average 0.7160000  0.7160000 0.8889723 0.8871116
## 
## Coverage of 95% prediction intervals:
##           EX.pred
## obs     0.9550562
## average 0.9206349

par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
     col=c("ID", "black", "grey"),
     ylim=c(-1,2),
     xlab="Observations", ylab="Predictions",
     main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
                                    xlab="Observations", ylab="Predictions",
                                    main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")

jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModB.jpeg"),
     width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
     col=c("ID", "black", "grey"),
     ylim=c(-1,2),
     xlab="Observations", ylab="Predictions",
     main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
                                      xlab="Observations", ylab="Predictions",
                                      main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()

## quartz_off_screen 
##                 2

3.2.5 Predictions for 2009-2020

What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites

x.B <- coef(est.denver.model.B, pars = "cov")$par
x.B

## [1] -8.693332 -7.074670 -7.568345 13.812075 -3.008518 -4.818215

E.B <- predict(denver.model.B, est.denver.model.B, LTA = T,
               transform="unbiased", pred.var=FALSE)
print(E.B)

## Prediction for STmodel.
## 
## Regression parameters:
##   0 Spatio-temporal covariate(s).
##   24 beta-fields regression parameters in x$pars.
## 
## Regression parameters are assumed to be known and
##  prediction variances do NOT include
##  uncertainties in regression parameters.
## 
## Prediction of beta-fields, (x$beta):
## List of 2
##  $ mu: num [1:69, 1:3] 0.339 0.182 0.183 0.225 0.19 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX: num [1:69, 1:3] 0.344 0.185 0.181 0.215 0.187 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Predictions for log-Gaussian field of type: unbiased 
## 
## Predictions for 624 times at 69 locations.
## List of 5
##  $ EX.mu     : num [1:624, 1:69] 1.45 1.48 1.6 1.53 1.55 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX.mu.beta: num [1:624, 1:69] 1.49 1.52 1.64 1.56 1.58 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX        : num [1:624, 1:69] 1.52 1.56 1.68 1.6 1.62 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX.pred   : num [1:624, 1:69] 1.53 1.56 1.69 1.61 1.62 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ log.EX    : num [1:624, 1:69] 0.396 0.419 0.493 0.445 0.456 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Variances have been computed.
## List of 2
##  $ log.VX     : num [1:624, 1:69] 0.05 0.0499 0.0498 0.0497 0.0496 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ log.VX.pred: num [1:624, 1:69] 0.0581 0.058 0.0579 0.0578 0.0577 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Mean squared prediciton errors NOT been computed.
## 
## 69 temporal averages have been compute.
## List of 1
##  $ LTA:'data.frame': 69 obs. of  4 variables:

week_preds <- as.data.frame(E.B$EX) %>%
  mutate(week = as.Date(rownames(E.B$EX)),
         year = year(as.Date(rownames(E.B$EX))))

week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]

box_preds <- week_preds %>%
  select(week, all_of(week_sites)) %>%
  #filter(week %in% week_weeks) %>%
  mutate(month = month(week),
         year = year(week)) %>%
  filter(year %in% c(2009:2020))

box_data <- box_preds %>%
  pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)

##       week                month             year        location        
##  Min.   :2009-01-05   Min.   : 1.000   Min.   :2009   Length:42364      
##  1st Qu.:2011-12-26   1st Qu.: 4.000   1st Qu.:2011   Class :character  
##  Median :2014-12-22   Median : 7.000   Median :2014   Mode  :character  
##  Mean   :2014-12-22   Mean   : 6.506   Mean   :2014                     
##  3rd Qu.:2017-12-18   3rd Qu.: 9.000   3rd Qu.:2017                     
##  Max.   :2020-12-07   Max.   :12.000   Max.   :2020                     
##                                                                         
##       pred       
##  Min.   :0.3674  
##  1st Qu.:1.1699  
##  Median :1.4733  
##  Mean   :1.5901  
##  3rd Qu.:2.0054  
##  Max.   :4.0348  
##  NA's   :1080

#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
  group_by(month) %>%
  summarize(median = median(pred, na.rm=T)) %>%
  arrange(desc(median))
box_summary

pred_box_plot <- ggplot(box_data) +
  geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
               show.legend = F) +
  #scale_color_viridis(name = "Month", discrete = T) +
  facet_wrap(. ~ year, scales = "free") +
  xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
  # theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot

4 Summary Table

Summary of model diagnostics

model	trend_data	cov_beta0	cov_beta1	cov_beta2	cov_nu_error	no_basis_fns	df_per_year	all_converge	cv_rmse_obs	cv_rmse_avg	cv_r2_obs	cv_r2_avg	cv_coverage_obs	cv_coverage_avg
Model A	BC + NO2 + PM2.5	iid	iid	NA	exp	1	4	TRUE	0.16	0.06	0.81	0.91	0.95	0.95
Model B	BC + NO2 + PM2.5	iid	iid	iid	exp	2	4	TRUE	0.15	0.06	0.82	0.89	0.96	0.92